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BACKGROUND 


In  developing  a  probabilistic  decompression  model  it  was  found  that  asymmetrical  gas 
kinetics  proved  to  fit  the  available  human  exposure  data  much  better  than  symmetrical 
exponential  kinetics  (1-3).  This  asymmetrical  model  used  exponential  kinetics  for  gas  uptake, 
but  when  tissue  gas  tension  exceeded  ambient  pressure  by  a  specified  amount  the  kinetics 
became  linear,  dramatically  slowing  the  rate  of  gas  elimination.  This  so-called  linear 
exponential  (LE)  model  was  initially  developed  for  computing  fixed  oxygen  partial  pressure 
decompression  tables  (4,5),  but  had  been  extended  to  air  in  some  exploratory  trials  (6). 

The  success  of  the  LE  model  led  us  to  seek  a  physiologic  rationale  that  might  explain 
its  performance.  One  rationale  was  that  the  change  in  kinetics  described  slowing  of  gas 
elimination  as  a  result  of  gas  phase  formation.  This  hypothesis  is  supported  qualitatively  by 
observations  that  the  slowing  of  inert  gas  elimination  during  decompression  (7-9).  Could  the 
slow  pace  of  gas  elimination  implied  by  the  success  of  the  LE  model  be  confirmed  by  a 
physiologically  and  physically  plausible  approximation  to  gas  elimination  in  the  presence  of 
bubbles? 

To  answer  this  question  we  wanted  to  model  gas  phase  growth  and  dissolution  in 
heterogeneous  architectures  representative  of  actual  tissue.  Previously,  a  Monte  Carlo 
simulation  of  diffusion  had  been  developed  in  a  tissue  model  to  explore  the  effect  of  tissue 
heterogeneity  and  microvascular  architecture  on  gas  exchange  (with  no  gas  phase  present) 
under  normobaric  conditions  (10-12).  As  a  first  step  towards  modeling  gas  phase  dynamics 
in  that  environment,  we  have  developed  a  model  of  bubble  evolution  during  decompression  in 
a  homogeneous  and  uniformly  perfused  tissue  by  extending  the  Monte  Carlo  methods  used  in 
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previous  simulations.  Using  a  homogenous  and  uniformly  perfused  tissue  allows  for 
comparisons  using  well  established  methods  of  modeling  bubble  evolution  (13).  However,  in 
principle  the  Monte  Carlo  approach  should  allow  modeling  in  architectures  of  any 
complexity.  This  report  presents  the  methods  employed  and  their  derivations,  results  of  tests 
examining  the  model's  ability  to  predict  bubble  radii  under  equilibrium  conditions,  a 
comparison  of  the  time  course  of  bubble  evolution  predicted  by  the  Monte  Carlo  model  with 
that  predicted  from  the  numerical  solution  of  a  partial  differential  equation  model,  and  the 
model's  performance  at  different  levels  of  precision.  Application  of  the  model  in  reports  to 
follow  will  help  us  see  if  gas  elimination  can  be  delayed  as  much  as  predicted  by  the  LE 
model. 

METHODS 

Overview 

This  section  provides  a  brief  description  of  the  bubble-liquid  module,  the  Monte  Carlo 
method,  initial  conditions,  and  the  simulation  structure.  Details  regarding  the  derivation  and 
implementation  of  the  techniques  described  are  presented  in  the  sections  that  follow. 

As  a  first  approximation  of  a  physiological  model  of  bubble  evolution  in  tissue,  we 
developed  a  model  of  a  single  spherical  bubble  evolving  in  a  homogenous  liquid  of  fixed  and 
finite  volume.  The  volume  of  the  liquid  surrounding  a  single  bubble  can  be  interpreted  to  be 
the  inverse  of  the  bubble  density  (number  of  bubbles  per  volume  of  liquid).  A  tissue  might 
consist  of  many  such  bubbles,  each  associated  with  its  own  liquid  volume  as  shown  in  Figure 
1.  The  system  is  open  to  mass  transfer. 
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FIGURE  1  -  Schematic  diagram  of  a  single  bubble-liquid  module.  The  shaded  central  region 
represents  the  bubble.  The  concentric  spheres  surrounding  the  bubble  are  the  shells  used  to 
simulate  the  concentration  gradient  in  the  liquid.  The  number  and  size  of  the  shells  and  the 
size  of  the  bubble  vary  according  to  the  parameters  of  the  model. 
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The  traditional  approach  to  modeling  the  temporal  and  spatial  evolution  of  inert  gas  in 
tissues  using  a  system  of  partial  differential  equations  becomes  intractable  when  the 
architecture  of  the  tissue  is  complex.  An  alternative  to  this  type  of  deterministic  modeling  is 
a  Monte  Carlo  simulation  (10-14).  The  essence  of  this  approach  involves  the  placement  of 
inert  gas  particles  in  a  bubble-liquid  module  and  allowing  them  to  take  random  steps  to 
simulate  diffusion  for  a  short  period  of  time.  At  the  end  of  the  time  period  the  distribution  of 
the  particles  is  used  to  calculate  the  number  of  moles  of  gas  in  the  bubble  and  liquid.  The 
new  bubble  volume  is  then  calculated  from  the  ideal  gas  law.  This  process  is  repeated  many 
times  to  obtain  the  time  course  of  bubble  evolution. 

We  start  our  simulation  with  the  liquid  saturated  with  inert  gas  at  one  ambient 
pressure  and  then  make  a  step  decrease  in  the  pressure.  Since  there  is  little  information 
about  the  initial  formation  of  bubbles  (15,16),  we  simply  assume  that  the  gas  phase  initially 
has  a  radius  close  to  the  critical  radius  computed  by  Ward  and  Tikuisis  (17,18).  We  then 
place  the  bubble  centered  in  the  liquid  region. 

We  obtain  the  value  of  the  critical  radius  from  calculations  of  equilibrium  radii  for  a 
closed  system  with  a  finite  gas  supply  carried  out  by  Ward  and  Tikuisis;  the  first  is  an 
unstable  radius,  referred  to  as  the  critical  radius  rc.  Gas  phases  smaller  than  rc  will  shrink 
and  disappear,  those  larger  will  grow  until  the  second  stable  equilibrium  radius,  re  is  reached. 
Bubbles  with  radii  larger  than  re  will  shrink  until  re  is  reached.  The  quartic  polynomial  from 
which  rc  and  re  are  obtained  is  given  in  Appendix  A.  The  specific  initial  condition  values  we 
used  in  our  system  test  are  given  in  the  section  entitled  "Model  Testing"  and  Table  1. 
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TABLE  1:  Model  Parameters  by  Equation  Number 


Eq  No.  Parameter 

Reference 

(1) 

Ptisoz  =  39.76  mmHg 

Ptiscoi  =  44-26  mmHg 

27 

(2) 

V,iq  =  0.000125  cm3,  0.0005  cm3  PDE 

Cliq  =  0.055140  moles/cm3 

Pjunho  =  1520  mmHg,  (2  ATA) 

Pvap  =  vapor  pressure  of  water  at  37  °C=  46.52  mmHg 
a  =  Ostwald  solubility  N2  =  0.0143 

Rg  =  62358  (cm3  mmHg)/(mole  K) 

T  =  310.1  K,  (37  °C) 

see  text 

1 /molecular  wt 

see  text 

27 

28 

36 

(3) 

pamb  =  760  mmHg,  (1  ATA) 

y  =  surface  tension  =  0.037515  mmHg-cm,  (50  dynes/cm) 
rbub  =  0.00065  cm,  0.00092  cm  PDE 

post-bubble  formation  boundary  condition  =  uniform  from  all 

see  text 

15 

see  text 

shells 

(23) 

DN2  =  0.5  X  10  5  cm2/sec,  0.75  X  10'5  cm2/sec  PDE 

see  text 

(26) 

grid  size  =  0.00025  cm 

[derived:  diffusion  step  size,  (a),  =  0.00043  cm] 

see  text 

(27) 

Hcycle  =  50 

[derived:  tcycle  =  0.3125  seconds] 

[derived:  tcycle  =  0.2100  seconds]  PDE 

see  text 

(57)  X 

=  2.5  minutes,  10  minutes  PDE 

see  text 

(62) 

Fi02  =  0.21,  (air) 
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The  amount  of  gas  required  to  fill  this  initial  bubble  is  taken  from  the  liquid  phase. 
This  can  be  done  either  by  removing  gas  uniformly  from  within  the  liquid,  leaving  no  initial 
gas  gradient  or  by  removing  gas  non-uniformly,  yielding  an  initial  gradient.  In  either  case 
we  keep  track  of  the  gas  concentration  with  a  series  of  concentric  shells  constructed  around 
the  bubble.  The  shells  expand  or  contract  with  the  expansion  or  contraction  of  the  bubble. 
The  fraction  of  total  gas  in  the  module  that  is  present  in  the  liquid  is  the  same  whether  or  not 
an  initial  gradient  is  present.  The  details  of  these  procedures  are  explained  in  the  sections 
entitled  "Calculation  of  Mole  Fraction  in  the  Bubble  and  Liquid  Shells"  and  "Simulation  of 
the  Gas  Gradient  in  the  Liquid." 

After  the  initial  placement  of  a  bubble  in  the  module  and  the  establishment  of  the  gas 
gradient  in  the  liquid,  the  simulation  begins  a  repeating  cycle  of  particle  placement,  random 
walks,  gas  wash-out,  and  adjustment  of  the  bubble  size  and  gas  gradient  in  the  liquid.  We 
refer  to  each  repetition  of  this  sequence  of  events  as  a  bubble  growth  cycle,  shown 
schematically  in  Figure  2.  The  simulation  consists  of  many  bubble  growth  cycles.  During 
each  cycle,  bubble  radius  is  kept  constant  and  no  gas  washout  occurs  while  the  gas  particles 
diffuse  throughout  the  module. 

Simulating  diffusion  using  a  Monte  Carlo  process  involves  many  particles, 
representing  idealized  gas  molecules,  being  placed  one  at  a  time  within  the  bubble-liquid 
module.  For  each  particle  the  location  of  the  placement  is  randomly  selected  to  be  in  either 
the  bubble  or  the  liquid,  based  on  the  fraction  of  total  gas  in  each  region.  The  derivation  of 
this  placement  procedure  is  contained  in  the  section  entitled  "Placement  of  Particles  in  the 
Bubble-liquid  Module"  and  Appendix  B. 


6 


Bubble  Growth  Cycle 


After  placement,  each  particle  walks  randomly  through  the  module  for  a  length  of 
time  specified  as  tcyde,  the  bubble  growth  cycle  time.  At  the  completion  of  the  random  walk, 
the  particle's  new  position  (i.e.,  which  shell  the  particle  is  in)  is  recorded.  Then  another 
particle  is  placed  and  walked  through  the  module.  This  is  repeated  for  all  of  the  particles  (at 
least  2  x  108  particles/(ml  of  liquid  in  the  module)).  The  details  of  this  procedure  are  in  the 
section  entitled  "Random  Walk  Procedure. " 

For  purposes  of  the  random  walk  procedure,  gas  particles  in  the  bubble  and  liquid 
phase  are  treated  differently.  If  the  gas  particles  are  in  the  bubble,  then  they  are  assumed  to 
be  uniformly  distributed  and  transition  from  the  bubble  to  the  liquid  can  only  occur  if  the 
particle  begins  its  step  from  within  a  shell  just  inside  the  bubble-liquid  boundary.  The 
probability  of  crossing  from  the  bubble  to  the  liquid  in  a  particular  step  is  the  product  of 
three  probabilities:  the  probability  of  being  in  the  transition  region,  the  probability  of  striking 
the  bubble-liquid  interface,  and  the  probability  of  crossing  the  interface  if  struck. 

If  in  the  liquid,  the  particles  are  allowed  to  cross  from  the  liquid  to  bubble  if  the 
bubble-liquid  boundary  is  struck.  If  a  particle  strikes  the  outer  liquid  boundary  during  a  step 
it  is  reflected  back  into  the  module.  This  simulates  a  zero  pressure  gradient  with  no  mass 
transfer  across  the  outer  tissue  boundary.  Details  are  in  the  section  entitled  "Derivation  of  the 
Transition  Probability." 

After  all  particles  have  been  placed  and  undergone  the  random  walk,  there  will  be  a 
new  particle  distribution.  The  number  of  particles  in  the  gas  phase  and  each  individual  shell 
is  divided  by  the  respective  volumes  of  the  gas  phase  or  shell,  to  get  the  mean  concentration 
in  the  gas  phase  and  in  each  shell.  Next,  the  concentration  in  each  shell  is  adjusted  to 
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FIGURE  3  -  Bubble  simulation  flow  chart.  rmin  is  equal  to  the  diffusion  step  size  or  critical 
radius,  whichever  is  larger. 
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-  Stepping  loop  flow  chart. 


This  algorithm  is  followed  whenever  the  particle  has 


a  chance  of  entering  the  bubble  during  a  bubble  growth  cycle  of  length  tcycle 


account  for  the  wash-in  or  wash-out  that  occurred  during  the  cycle  time.  Thus,  the  gas 
concentration  gradient  in  the  liquid  is  established  for  the  beginning  of  the  next  cycle.  The 
bubble  radius  is  adjusted  according  to  the  number  of  moles  in  the  gas  phase  assuming  the 
ideal  gas  law  as  explained  in  the  section  entitled  "Calculation  of  the  New  Bubble  Radius." 
Once  this  new  radius  is  established  the  thickness  of  the  shells  is  adjusted.  The  new  shell 
radii  and  the  mean  gas  concentration  in  each  of  the  shells  are  the  initial  conditions  for  the 
next  cycle.  Details  of  this  process  are  given  in  the  section  entitled  "Simulation  of  the  Gas 
Gradient  in  the  Liquid." 

In  summary,  we  begin  with  a  liquid  saturated  with  inert  gas  at  some  ambient  pressure, 
reduce  the  ambient  pressure,  and  insert  a  spherical  bubble  with  a  radius  just  above  the  critical 
radius  for  the  system.  For  a  short  time  interval  we  fix  the  bubble  size  while  gas  in  the 
bubble-liquid  module  is  redistributed  using  a  Monte  Carlo  simulation  of  diffusion.  At  the  end 
of  this  time  interval,  the  bubble  size  and  gas  gradient  in  the  liquid  are  adjusted  based  on  this 
new  particle  distribution.  The  next  cycle  then  begins,  and  this  process  continues  indefinitely. 
Figures  3  and  4  present  flow  charts  that  summarize  the  model's  logical  structure.  Details  of 
the  derivation  and  implementation  of  the  method  follow. 

Calculation  of  Mole  Fraction  in  the  Bubble  and  Liquid  Shells 

The  number  of  moles  of  inert  gas  dissolved  in  the  liquid  at  the  start  of  the  simulation 
before  a  bubble  forms  (n^J  is  calculated  from  the  saturated  volume  of  inert  gas  at  the  ambient 
pressure  just  before  the  pressure  reduction.  In  living  organisms  the  metabolic  gases  oxygen 
and  carbon  dioxide  are  present  and  must  be  taken  into  account.  We  first  define  the  partial 
pressure  from  metabolic  gases  (Pmg)  to  be, 
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(1) 


mg 


=  iV 


+P. 


°2 


tlSr 


L0i 


where  Ptis02  =  partial  pressure  of  tissue  oxygen  (mmHg) 

PtisC02  =  partial  pressure  of  tissue  carbon  dioxide  (mmHg) 

For  simulations  not  involving  metabolic  gases  pmg  is  assumed  to  be  0.  i\ot  is  then  calculated 

as, 


V  C  P 

_  "  tiq  liq  rljq 

"'tot 


Kt 


H 


(2) 


where, 

Vjiq  =  volume  of  liquid,  (cm3) 

Cliq  =  molar  concentration  of  the  solvent,  (moles/cm3) 

Pliq  =  pressure  of  inert  gas  in  the  liquid,  (mmHg) 

(^ambO'^vap’^mg) 

Pambo  =  ambient  pressure  before  pressure  reduction,  (mmHg) 

Pvap  =  vapor  pressure  of  water,  (mmHg)  (constant) 

Kh  =  Henry's  coefficient  =  PH/([(Rg-T/(VL-a-PH))  +  l]'1),  (mmHg)(35) 

VL  =  molar  volume  of  solvent  (1/Cliq),  (cm3/mole) 
a  =  Ostwald  solubility  coefficient 
PH  =  760  mmHg 

Rg  =  Universal  gas  constant  in  (cm3  mmHg)/(mole  K) 

T  =  temperature  in  Kelvin 

Once  a  bubble  is  placed  in  the  liquid,  the  liquid  is  assumed  to  remain  at  constant 
volume  as  the  bubble  expands.  All  the  gas  in  the  bubble  comes  from  the  liquid  phase, 
keeping  the  total  amount  of  gas  in  the  module  constant  during  a  bubble  growth  cycle. 
Assuming  an  ideal  gas  and  negligible  liquid  elastic  forces,  the  number  of  moles  of  inert  gas 
in  the  bubble  (n^)  is  determined  by, 
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where, 

Pbub  =  inert  gas  pressure  in  the  bubble  =  Pamb-Pvap-Pmg+Pr,  (mmHg) 

Pam5  =  ambient  pressure  after  pressure  reduction,  (mmHg) 

Py  =  surface  tension  pressure  =  (2*y)/rbub,  (mmHg) 
y  =  surface  tension,  (mmHg-cm) 

[y  (mmHg-cm)  =  y  (dynes/cm)  /  1332.8  dynes/(mmHg-cm2)] 
rbub=  bubble  radius,  (cm) 

Vbub  =  bubble  volume,  (cm3).  (Vbub  is  the  volume  of  a  bubble  of  radius  rbub  selected  to  be 
just  above  the  critical  radius  for  a  pressure  reduction  from  Pamb0  to  Pamb.) 

Because  placement  of  particles  in  the  module  for  the  Monte  Carlo  simulation  of 

diffusion  depends  on  the  relative  amount  of  gas  in  each  region  of  the  bubble-liquid  module,  it 

is  necessary  to  calculate  the  mole  fraction  of  gas  in  each  region.  The  initial 

mole  fraction  in  the  bubble  is  simply. 


Xbub0 


n 


bub 


n, 


tot 


(4) 


and  the  initial  mole  fraction  in  the  liquid  is, 


n 


bub 


n. 


tot 


The  initial  mole  fraction  in  each  shell,  xshelli0,  is  determined  from  the  equation, 


(5) 
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At  the  start  of  the  simulation,  the  number  of  moles  of  gas  in  each  shell,  n^,  is  determined 
by  the  selection  of  the  gas  gradient  in  the  liquid.  We  have  provided  for  one  of  two 
alternative  gas  distributions:  the  gas  for  the  bubble  can  be  taken  uniformly  from  all  of  the 
shells  in  the  liquid  so  that  there  is  no  initial  gradient  in  the  liquid,  or  the  inner-most  shells 
can  be  depleted  of  the  gas  that  makes  up  the  bubble  so  that  partial  pressure  equilibrium 
across  the  boundary  is  maintained.  (These  two  choices  represent  the  largest  and  smallest 
gradients  for  gas  diffusion.  We  may  then  evaluate  the  effect  of  these  two  extreme  post¬ 
bubble  formation  boundary  conditions  on  bubble  evolution.)  In  the  uniform  distribution  case, 
the  number  of  moles  in  each  shell  is  the  same  and  is  simply, 


shell‘  number  of  shells 


In  the  other  case,  the  procedure  for  calculating  the  number  of  moles  in  each  shell  is  more 
complex.  Gas  is  taken  from  the  innermost  shell  until  the  partial  pressure  is  equal  to  the 
calculated  partial  pressure  for  inert  gas  in  the  bubble.  If  additional  gas  is  required  to  fill  the 
bubble,  the  same  procedure  is  applied  to  succeeding  shells  until  the  total  gas  needed  has  been 
taken.  Appendix  D4  contains  the  FORTRAN  code  that  implements  this  procedure. 

After  the  first  bubble  growth  cycle,  the  particles  will  have  a  different  distribution  than 
before  their  random  walks  simulating  diffusion  began.  This  new  distribution  forms  the  basis 
for  placement  to  begin  the  next  bubble  growth  cycle.  Consequently,  it  is  necessary  to 
calculate  the  mole  fraction  in  each  region  of  the  module  at  the  end  of  the  cycle.  In  this  case, 
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the  mole  fraction  in  each  region  is  determined  by  the  fraction  of  particles  that  are  in  the 
region  at  the  end  of  a  particular  bubble  growth  cycle,  which  is  given  by, 


Xbub 


number  of  particles  in  the  bubble 
total  number  of  particles  in  the  module 


(8) 


The  mole  fraction  in  the  liquid  is, 


number  of  particles  in  the  liquid  _t 

Xu  '  1  X huh 

q  total  number  of  particles  in  the  module 


(9) 


and  the  mole  fraction  in  each  shell  is, 


number  of  particles  in  shell { 
shell‘  total  number  of  particles  in  the  module 

Placement  of  Particles  in  the  Bubble-liquid  Module 

In  order  to  obtain  the  correct  distribution  of  particles  after  they  are  subject  to  diffusion 
and  wash-out,  the  placement  of  particles  at  the  beginning  of  each  bubble  growth  cycle  cannot 
be  done  haphazardly,  but  must  follow  certain  rules. 

First,  a  random  number  between  0  and  1  is  compared  with  the  mole  fraction  in  the 
bubble.  If  the  random  number  is  less  than  the  mole  fraction  the  particle  is  assigned  to  the 
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bubble,  otherwise  it  is  assigned  to  the  liquid.  The  liquid  is  divided  into  concentric  shells  of 
equal  width,  which  are  used  to  simulate  the  concentration  gradient  of  inert  gas  in  the  liquid. 
Within  the  liquid,  the  probability  of  placement  in  each  shell  is  proportional  to  the  mole 
fraction  of  gas  in  the  shell.  In  order  to  weight  the  particle  placement  by  the  mole  fraction  in 
each  shell,  a  random  number  between  0  and  1  is  compared  with  the  serially  accumulating 
mole  fraction  beginning  from  the  inner  most  shell.  At  the  point  the  random  number  is  less 
than  that  value,  that  shell  is  selected  for  particle  placement.  The  smooth  gradient  in  the 
liquid  is  approximated  by  the  series  of  shells  with  different  mean  concentrations  in  each  shell, 
forming  a  step  gradient.  However,  the  concentration  is  uniform  within  each  shell,  implying 
that  the  probability  of  a  particle  being  placed  in  any  two  regions  of  equal  volume  within  a 
given  shell  must  be  the  same.  Within  a  shell,  our  spherical  coordinate  system  requires  that 
more  particles  be  placed  farther  from  the  bubble  because  a  spherical  sub-shell  far  from  the 
bubble  will  occupy  more  volume  than  one  of  equal  width  nearer  the  bubble. 

We  use  the  spherical  coordinates  p,  0,  and  (j)  to  place  the  particles  within  a  shell, 
where  p  is  the  radial  coordinate,  0  the  azimuthal  angle  between  0  and  7t,  and  (j)  orthogonal  to 
0,  between  0  and  271.  In  order  to  ensure  the  proper  values  for  these  variables,  we  must 
transform  the  uniformly  distributed  random  numbers  (u)  between  0  and  1  generated  by  the 
computer.  We  do  this  using  the  transformations  hp,  he,  and  h^  derived  in  Appendix  B  (19). 


p  =  hp(u )  =  Ru 


(ID 


where  R  is  the  radius  of  the  sphere.  However,  since  placement  is  to  occur  in  a  shell  and  not 
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an  entire  sphere,  pmust  be  adjusted  accordingly, 


,  3  ,3  3  s>.  ' 

P  =  (r{_i  +  uirt  -rM))- 


(12) 


with  p  between  the  outer  (r;)  and  inner  (r^)  shell  radii. 


0  =  hJu )  =  arccos(l  -2u) 


(13) 


and 


<j)  =  h^(u)  -  2k -u 


(14) 


These  functions  map  the  uniformly  distributed  random  variable  u  between  0  and  1  into  the 
values  p,  0  and  (|)  which  are  uniformly  distributed  by  volume  within  a  shell.  These 
coordinates  are  translated  to  cartesian  coordinates, 


x  =  p-sin(0)-cos((J)) 


(15) 


y  -  p-sin(0)sin(<J)) 


(16) 


z  =  p-cos(0) 


(17) 


Random-walk  Procedure 

The  simplest  way  of  simulating  diffusion  with  a  Monte-Carlo  process  is  to  translate 
from  the  spherical  coordinate  system  of  Figure  1  into  a  cartesian  coordinate  system  using 
Equations  15-17.  A  random  number  generator  would  be  used  to  specify  positive  or  negative 
movement  in  the  X,  Y,  and  Z  directions  (10-12),  with  the  length  of  the  movement  in  each 
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dimension  being  equal  to  a  cartesian  grid  unit.  The  particle  would  end  up  on  one  of  the  8 
corners  of  the  cube  with  sides  2  grid  units  in  length  centered  on  the  original  position.  Thus 
the  effective  distance  traveled,  which  we  refer  to  as  the  diffusion  step  size  (a),  would  be 
V 3 -(length  of  the  grid  unit).  This  process  would  be  repeated  for  a  predetermined  number  of 
steps  and  would  result  in  the  generation  of  the  particle's  path  through  the  module. 

Unfortunately,  this  method  requires  large  amounts  of  computer  execution  time.  We 
took  advantage  of  the  fact  that  the  region  outside  the  bubble  is  isotropic  to  diffusion  by 
calculating  the  distribution  of  distances  a  particle  will  travel  in  a  relatively  large  amount  of 
time  under  the  influence  of  diffusion  alone.  We  refer  to  this  as  the  Bigstep  distribution.  This 
distribution  is  computed  as  follows.  If  an  amount  of  substance  M  is  deposited  at  a  point  in 
an  infinite  volume,  then  the  concentration  C  at  a  distance  r  from  the  point  source  under  the 
influence  of  diffusion  alone  is  given  by  Crank  (20)  as, 

-r2 

C-  M  (18) 

(4ti  Dif2 


By  dividing  both  sides  of  this  equation  by  M  we  obtain  the  probability  density  function  (PDF) 
of  the  distance  r,  a  particle  will  travel  in  time  t  from  a  point  in  the  volume, 


Ar,t)  = 


1 


(4; iDt)312 


•e 


-r 
4 Dt 


(19) 


The  cumulative  distribution  function  (CDF)  can  then  be  obtained  by  integrating  the  PDF  over 
the  volume, 
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(20) 


r  2«  n  -r2 


F(r,t) = - f  f  fe  4Dtr2-sin0ddd<t>dr 


which  reduces  to, 


F(r,r)= 


471 

(4tiD03/2 
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fr2-e  4Dtdr 
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(21) 


r 

z= - 

If  we  let  V/25t  this  becomes, 


F(z)- 


N 


it  i 


z2-e  2  dz 


(22) 


Since  a  closed  form  solution  of  this  integral  cannot  be  obtained,  we  integrate  it  numerically  to 
obtain  F(z),  which  is  the  probability  a  particle  will  travel  the  normalized  distance  z  in  time  t. 


We  set  the  upper  bound  of  the  CDF  integral  at  5  normalized  distance  units  ( 


572 Dt 


) 


because  99.99 %  of  the  particles  will  travel  less  than  or  equal  to  this  distance  in  time  t. 
Graphs  of  the  PDF  and  numerically  obtained  CDF  are  shown  in  Figure  5. 

By  applying  the  inverse  transformation  method  (13)  to  the  CDF  by  randomly 
generating  a  number  between  0  and  1,  we  can  select  the  normalized  distance  z  a  particle 
travels.  We  do  this  by  discretizing  F(z)  into  increments  of  0.0001  and  assign  to  each  the 
value  of  z  needed  to  obtain  it.  These  z  values  are  stored  in  a  10,000-element  array.  We 
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Normalized  Distance 


FIGURE  5  -  Probability  and  cumulative  distribution  functions  of  normalized  diffusion 
distances  in  three  dimensions.  The  solid  line  represents  the  PDF  and  the  dashed  line 
represents  the  CDF. 
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then  randomly  generate  a  value  for  F(z)  and  go  to  the  array  to  find  the  associated  value  of  z. 


r 

z= - 

Since  ,  we  must  either  select  a  value  for  r  and  compute  t,  or  conversely  specify 

t  and  compute  r.  We  elect  to  keep  the  diffusion  time  increment,  t,  constant  and  compute  the 
distance  traveled  as, 

r  -  zfiDt  <23> 

The  value  we  choose  for  t  is  limited  by  the  architecture  of  the  module.  One 
constraint  placed  on  the  procedure  is  that  we  do  not  want  the  particle  to  cross  the  liquid-gas 
boundary  during  a  bigstep.  Once  a  particle  enters  the  bubble  it  must  be  given  the  opportunity 
to  exit,  diffuse  in  the  liquid  and  reenter  the  bubble,  perhaps  on  several  occasions.  If  we 
allowed  the  particle  to  enter  the  bubble  on  a  bigstep  this  opportunity  would  not  be  provided 
and  we  would  end  up  with  the  wrong  distribution  of  particles  between  the  phases  of  the 
module. 

We  apply  this  constraint  by  first  requiring  that  the  bigstep  only  be  allowed  for 
particles  outside  of  the  first  shell.  Secondly,  we  require  that  the  probability  of  the  particle 
reaching  the  bubble  boundary  during  a  bigstep  be  very  small  (<  0.0001).  Since  the  upper 
limit  of  the  integral  in  Equation  22  is  5  normalized  distance  units,  the  probability  of  attaining 
a  distance  of  greater  than  or  equal  to  5  normalized  distance  units  is  less  than  0.0001.  If  the 

r0 

particle  is  a  distance  r0  from  the  boundary,  then  as  long  as  - ^5  the  particle  has  a 

yJWt 
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very  small  chance  (<  0.0001)  of  reaching  the  bubble  boundary  during  time  t. 

We  define  the  largest  distance  a  particle  can  travel  in  a  straight  line  on  a  single 
diffusion  step  to  be  the  mean  squared  distance.  The  mean  squared  distance  (MSD)  a  particle 
travels  under  the  influence  of  diffusion  in  three  dimensions  is  (20), 

MSD  =  6Dt  (24> 

where  D  is  the  diffusion  coefficient  and  t  is  the  time  allowed  for  diffusion.  At  the  beginning 
of  the  simulation,  we  assign  a  value  for  the  mean  squared  distance  on  a  single  diffusion  step 
to  be  a2  where, 

a2  =  6Dt  (25) 

step 


We  solve  Equation  25  for  t^,  to  get  the  time  to  travel  the  distance  (a)  that  corresponds  to  a 
single  diffusion  step. 


step 


(26) 


We  also  assign  a  value  to  the  number  of  single  steps  each  particle  will  undergo  during  each 
bubble  growth  cycle,  n^,  and  compute  the  time  for  each  bubble  growth  cycle,  tcycle,  as 


cycle 


= 

cycle  step 


(27) 


The  time  tcycle  is  the  longest  time  allowed  for  a  bigstep,  so  if 


2:5  then  there  is 
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almost  no  chance  the  particle  will  reach  the  boundary  during  time  tcyc]e.  In  this  case  the 
distance  rx  actually  traveled  during  the  time  tcycle  can  be  calculated  directly  from  Equation  23 

aS  rT  =  cycle' 

The  procedure  can  be  illustrated  with  an  example.  If  the  computer  generated  random 
number  is  0.12145673,  we  multiply  it  by  10,000  and  round  to  obtain  the  array  index,  1,215. 
The  associated  z  value  of  0.823  (obtained  from  a  table  look-up)  combined  with  D  of  0.5X1  O'5 
cm2/sec  and  t  equal  to  tcycle  of  0.3125  sec  in  Equation  23,  gives  the  actual  distance  traveled  of 
0.00173  cm  (17.3  microns). 
r0 

If  <5  then  the  time  allowed  for  bigstep  is  reduced  as  follows.  First,  the 

cycle 


number  of  single  steps  of  length  a  in  a  direct  path  to  the  boundary  is  calculated  as, 


Hstep  =  a 


(28) 


then 


*  _  y. 

lbigstep  ~  rlstep  lstep 


(29) 


and 


rT  Zyj2Dt bigstep 


(30) 
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The  probability  of  all  of  the  random  walk  steps  occurring  in  the  same  direction  is  small,  so 
by  computing  the  distance  traveled  in  time  tbigstep  using  the  Bigstep  distribution,  the  particle 
will  have  almost  no  chance  of  stepping  distance  r0  and  entering  the  bubble. 

Although  this  is  the  method  we  have  used  in  the  current  version  of  our  program,  a 
more  consistent  approach  would  be  to  calculate  t^p  from  Equation  23  after  transforming  the 
radial  distance  to  the  bubble  boundary  into  5  normalized  distance  units. 


1 bigstep 


(31) 


Then  by  randomly  generating  z  as  previously  described,  and  using  tbigstep  from  Equation  31, 
the  distance  traveled  can  be  calculated  using  Equation  30. 

This  approach  guarantees  that  the  particle  will  have  less  than  a  0.01  %  chance  of 
entering  the  bubble  on  a  bigstep. 

The  direction  of  the  step  is  obtained  by  randomly  generating  two  orthogonal  angles,  0 
and  fusing  Equations  13  and  14.  Thus,  a  particle  will  be  placed  randomly  on  the  surface  of 
a  sphere  centered  at  its  initial  location  with  a  radius  obtained  randomly  from  the  Bigstep 
distribution. 

For  those  particles  with  tbjgstep<tcycle,  the  particles  must  be  given  the  opportunity  to 
take  individual  diffusion  steps  as  described  in  the  beginning  of  this  section  to  complete  the 
bubble  growth  cycle.  The  time  remaining  in  the  bubble  growth  cycle  after  taking  the  bigstep 
is  obtained  and  the  number  of  single  diffusion  steps  is  calculated  as, 
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t cycle  ^ bigstep 


(32) 


aicy  , 

Lstep 

If  the  particle  has  crossed  the  bubble  boundary  and  entered  the  gas  phase,  a  different 
stepping  procedure  is  invoked.  We  can  calculate  the  probability  the  particle  will  be  available 
to  exit  from  the  bubble  during  a  single  step,  which  is  equal  to  the  product  of  three  individual 
probabilities:  the  probability  of  being  in  the  transition  region,  the  probability  of  striking  the 
bubble-liquid  interface  from  the  transition  region,  and  the  probability  of  crossing  the  interface 
if  struck.  If  a  random  number  generated  by  the  computer  is  less  than  the  product  of  these 
probabilities,  the  particle  is  placed  in  the  transition  region  uniformly  by  volume  as  described 
in  the  section  on  placement  of  particles.  It  is  then  allowed  to  take  a  single  step  to  the  surface 
of  a  sphere  of  radius  (a),  centered  at  its  location  in  the  transition  region.  If  this  step  does  not 
take  the  particle  outside  the  boundary,  it  is  "returned"  to  the  well-mixed  bubble  region  and 
given  another  opportunity  to  be  available  for  exiting.  This  cycle  is  repeated  until  the  particle 
steps  out  of  the  bubble  or  the  number  of  steps  allowed  in  a  bubble  growth  cycle  is  reached. 
The  derivation  of  this  procedure  is  described  in  the  section  "Derivation  of  the  Bubble 
Transition  Probability." 

Any  particle  that  reaches  the  outer  boundary  of  the  liquid  module  is  reflected  inward 
to  simulate  a  particle  entering  from  an  adjacent  bubble-liquid  module.  A  slight  error  is 
introduced  because  packed  spheres  are  not  space  filling,  and  particles  can  never  enter  this 
intersphere  region. 

Simulation  of  the  Gas  Gradient  in  the  Liquid 

In  order  to  simulate  a  concentration  gradient  within  the  liquid,  "shells"  consisting  of 
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concentric  spheres  are  placed  around  the  bubble.  At  the  start  of  each  bubble  growth  cycle 
the  width  of  each  shell  is  k^lDt  k  , where  D  is  the  diffusion  coefficient,  tcycle  is  cycle  time, 

and  k  is  a  constant  between  0  and  1 .  The  root  mean  square  distance  the  particle  will  travel 
in  tcycle  is  y/6 Dt cycle-  Since  a  particle  will  travel  this  average  distance  in  tcycle,  the  assumption 

that  the  shells  are  well  mixed  is  consistent  with  a  maximum  shell  width  of  j2Dt ^  ■  This 

choice  for  shell  width  is  also  convenient,  since  it  is  proportional  to  the  normalized  distance 
unit  of  the  Bigstep  distribution.  For  greater  resolution  of  the  gradient,  k  can  be  selected 
smaller  than  1. 

After  the  random  walk  portion  of  the  growth  cycle  the  bubble  will  change  size.  This 
will  result  in  a  change  in  the  thickness  of  the  shells  surrounding  the  bubble.  For  example,  if 
the  bubble  grows  the  inner  and  outer  radii  of  each  shell  increase,  which  makes  the  shell 
thickness  become  smaller  since  the  shell  volume  is  held  constant  during  expansion.  For  a 
growing  bubble  the  shell  width  would  eventually  become  very  small.  In  order  to  adjust  for 

this,  a  new  set  of  shells  of  width  k-^jlDt^  is  created  at  the  beginning  of  each  cycle.  These 

new  shells  will  overlap  the  expanded  or  contracted  shells  of  the  previous  cycle  as  shown  in 
Figure  6.  While  they  are  the  same  thickness  as  the  previous  set  of  shells  before  the  bubble 
changed  size,  they  will  not  have  the  same  volume  as  the  expanded  or  contracted  shells  of  the 
previous  cycle.  Particles  are  placed  using  the  expanded/contracted  shells  of  the  previous 
cycle.  The  new  set  of  shells  is  used  to  record  the  final  position  of  the  particles  at  the  end  of 
their  random  walks. 
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FIGURE  6  -  Schematic  diagram  of  cross-sectional  view  of  changes  in  shells  during  bubble 
expansion.  In  the  first  frame,  the  dashed  lines  represent  the  location  of  the  shell  boundaries 
prior  to  expansion.  After  the  bubble  expands,  the  shell  boundaries  also  expand  to  the  solid 
lines  so  that  shell  volumes  are  the  same.  In  the  second  frame,  a  new  set  set  of  shells,  all 
with  the  same  width,  are  overlayed  on  the  expanded  shells  of  frame  1 .  Placement  of 
particles  in  these  new  shells  is  weighted  by  the  expanded  shells. 
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Derivation  of  the  Bubble  Transition  Probability 

In  our  simulation  the  gas  phase  is  assumed  homogeneous  and  a  particle  in  the  gas 
phase  has  an  equal  probability  of  being  anywhere  in  the  bubble.  If  a  particle  initially  in  the 
liquid  phase  encounters  the  bubble  boundary,  it  is  allowed  to  cross  the  boundary.  Since  a 
single  diffusion  step  is  the  largest  straight  line  distance  a  particle  is  allowed  to  travel 
(representing  the  Mean  Squared  Distance),  a  particle  is  not  allowed  to  enter  and  again  exit 
the  bubble  during  a  single  diffusion  step. 

If  a  particle  is  placed  in  the  bubble,  we  then  compute  a  probability  that  it  will  exit  on 
a  particular  step.  To  compute  this  probability  we  first  define  an  inner  transition  region  within 
the  bubble  which  is  a  spherical  shell  of  thickness  {a),  whose  outer  boundary  is  the  gas-liquid 
interface.  Only  particles  placed  in  this  inner  transition  region  are  allowed  to  exit  the  bubble 
(with  a  certain  probability),  particles  not  in  the  transition  region  cannot  exit  the  bubble. 

The  probability  of  exiting  the  bubble  is  the  product  of  three  individual  probabilities: 
the  probability  of  being  in  the  inner  transition  region,  the  probability  of  striking  the  bubble- 
liquid  interface  from  the  transition  region  and  the  probability  of  crossing  the  interface  if 
struck. 

The  probability  of  being  in  the  inner  transition  region  is  simply  the  ratio  of  the  inner 
transition  region  volume  to  the  bubble  volume, 


P 


trans 


Ri-jR-a? 
R 3 


(33) 


Once  in  the  inner  transition  region  the  probability  of  striking  the  interface  is  derived 
as  follows.  The  particle  can  step  from  point  P  in  the  inner  transition  region  to  any  point  a 
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FIGURE  7  -  Schematic  diagram  of  transition  region  inside  the  bubble.  The  upper  part  of  this 
illustration  represents  a  cross  section  through  the  bubble  centered  at  point  C  with  radius  R. 
The  transition  region  (T.R.)  extends  a  distance  a  inward  from  the  bubble  boundary.  The 
lower  part  of  the  illustration  is  a  cross  section  through  the  transition  sphere  of  radius  a 
centered  at  point  P  in  the  transition  region.  The  point  I  is  as  the  intersection  of  the  transition 
sphere  and  the  bubble  cross  section,  r  is  the  distance  of  point  P  from  the  surface  of  the 
bubble. 
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distance  (a)  away,  as  illustrated  in  Figure  7.  The  probability  of  striking  the  interface  (Pint)  is 
the  fractional  surface  area  of  the  sphere  of  radius  (a)  that  lies  outside  of  the  bubble.  The 
fraction  of  this  small  sphere  outside  of  the  bubble  is  given  by, 

P  (e)=(1-cos(9)>  04) 

mt  2 


where  0  is  the  angle  between  line  CP  and  line  PI,  and  I  is  a  point  of  intersection  of  the  two 
spheres.  Using  the  law  of  cosines,  we  obtain  the  following  expression  for  cos(0), 


cos(0)= 


2  rR-r2-a2 
2 a(R-r) 


(35) 


where  r  is  the  distance  of  point  P  from  the  surface  of  the  bubble.  The  probability  that  a 
particle  strikes  the  bubble  interface  from  point  P  is  then. 


PJr) 


r2  -2(a+R)r+a2  +2aR 
4a(i?-r) 


(36) 


To  get  the  probability  of  striking  the  interface  from  anywhere  in  the  transition  region,  we 
need  to  integrate  this  expression  over  the  entire  region,  weighted  appropriately. 


P.  = 


mt 


fPJrX4n(R-rf)dr 
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f  (4n(R-r)2)dr 


(37) 


Substituting  in  the  expression  for  Pint(r),  and  integrating  yields 
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P  =  a(12/g2  a2)  (38) 

M  16(/?3-(/?-a)3) 

If  the  particle  strikes  the  interface,  the  probability  of  crossing  the  interface  is  simply  the 
Ostwald  solubility  coefficient, 

p  =  a  (39) 

r cross  v  7 


The  individual  probabilities  can  be  combined  to  obtain  the  probability  of  exiting  the 
bubble  on  a  single  diffusion  step  pexit, 

Pexit  =  P cross'P trans’P int  (4( 

which  reduces  to, 


a(12R2-a2)  (dU 

Peat  =  a~ - ; — ~  <41> 

16/?3 

Verification  of  the  Transition  Algorithm  -  Because  the  transition  algorithm  extends  methods 
originally  developed  assuming  diffusivity  was  approximately  the  same  in  high  and  low 
solubility  regions,  we  felt  obligated  to  present  an  alternative  argument  in  support  of  using  this 
method  in  our  simulation.  The  second  approach  to  calculating  the  probability  of  leaving  a 
gas  bubble  begins  by  defining  several  probabilities. 

We  define  an  outer  transition  region,  also  of  thickness  (a),  in  the  liquid  immediately 
adjacent  to  the  bubble  as  shown  in  Figure  8.  P(t  |g)  is  the  probability  of  entering  the  outer 
transition  region  given  that  the  particle  is  in  the  gas  phase.  P(g  |t)  is  the  probability  of 
entering  the  gas  phase  given  that  the  particle  is  in  the  outer  transition  region.  If  P(t)  is  the 
probability  of  being  in  the  outer  transition  region,  and  P(g)  is  the  probability  of  being  in  the 
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FIGURE  8  -  Schematic  diagram  of  transition  region  outside  the  bubble.  The  upper  part  of 
this  illustration  represents  a  cross  section  through  the  bubble  centered  at  point  C  with  radius 
R.  The  transition  region  (T.R.)  extends  a  distance  a  outward  from  the  bubble  boundary. 

The  lower  part  of  the  illustration  is  a  cross  section  through  the  transition  sphere  of  radius  a 
centered  at  point  P  in  the  transition  region.  The  point  I  is  at  the  intersection  of  the  transition 
sphere  and  the  bubble  cross  section,  z  is  the  distance  of  point  P  from  the  surface  of  the 
bubble. 
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gas  phase  then,  at  equilibrium,  particles  must  be  transferring  from  liquid  to  gas  at  the  same 
rate  as  from  gas  to  liquid,  so  we  have, 


P(t\g)P(g)  =  P(g\t)P(t)  (42) 

Also,  at  equilibrium,  since  the  partial  pressures  of  the  gas  must  be  equal  on  either  side,  the 
probabilities,  P(i)  of  being  in  phase  i  with  volume  Vj  are  related  to  the  Ostwald  solubility 
coefficient  a  by, 


P(g)  _  vs 


(43) 


Pit)  Vta 

These  two  equations  may  be  solved  for  P(t  |g)  as  a  function  of  P(g  |t),  the  volumes, 
and  the  solubility  coefficient. 

P(g  |t)  is  a  geometric  factor  that  we  calculate  as  follows.  Imagine  a  particle  within  the 
outer  transition  region,  a  shell  of  thickness  (a)  surrounding  the  bubble  of  radius  R,  and 
located  a  distance  z  <  a  from  the  gas  phase  boundary  as  in  Figure  7.  The  next  diffusion  step 
of  size  (a)  defines  a  sphere  of  radius  (a).  The  probability  that  the  particle  will  enter  the  gas 
phase  is  proportional  to  the  surface  area  of  this  small  sphere,  which  lies  within  the  bubble.  If 
the  sphere  of  radius  (a)  intersects  the  bubble  at  an  angle  0  from  the  line  normal  to  the  surface 
of  both  spheres,  then  0  satisfies, 


Cos(0)  = 


az+2Rz+z: 
2 a(R+z) 


(44) 


and  P(g  1 1)  is  given  by  the  fraction 
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Pfelo  -  ^ 


(45) 


Then  P(t  |g)  can  be  obtained  in  terms  of  a,  R,  z,  and  a  by  substitution  into  the  equation  for 
P(t  |g)  above  with  the  expression  for  P(g  |t)  and  writing  Vg  in  terms  of  R. 

P(t\g)  =  cc(l-cos6)((R+a)3-R3)  (46) 

2  R3 

Using  the  expression  for  cos0,  we  obtain, 

p(,i  )  =  a[2a(R+z)-(a2+2Rz+z2)][(.R+a)3-R 3]  (47) 

4  aR\R+z) 

Notice  that  as  z  tends  to  (a),  P(t  |g)  tends  to  zero.  We  interpret  this  to  mean  that  exits 
into  the  transition  region  do  not  occur  with  a  frequency  representing  the  volume,  but  are 
biased  in  closer  to  the  bubble.  The  cumulative  probability  of  exit  is. 


faP(t\g)4n(R+z)2dz 
F(t\g)  =  - 

f  4n(R+z)zdz 
Jo 


(48) 


which  integrates  directly  to  give, 


F(t\g)  =  a 


a(12R2-a2) 
1 6R3 


(49) 


This  is  equivalent  to  the  expression  Pexit  used  in  the  transition  algorithm. 

Calculation  of  the  New  Bubble  Radius 

After  completing  a  bubble  growth  cycle,  the  location  of  all  particles  placed  during  the 
cycle  is  known.  The  mole  fraction  in  the  bubble,  xbub  is  taken  to  be  the  ratio  of  the  number 
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of  particles  ending  in  the  bubble  and  the  total  number  of  particles  placed  during  the  cycle 
from  Equation  8.  The  number  of  moles  of  gas  in  the  bubble  is  then  calculated  as, 


nbub  Xbubntot 


(50) 


Assuming  the  ideal  gas  law  and  using  the  expression  for  bubble  pressure  derived  in  equation 
3  gives, 


' bub 


(51) 


This  reduces  to  the  cubic  polynomial, 


YrL.-»»«r  -  0  (52) 

The  real  positive  root  of  this  polynomial  is  the  new  bubble  radius.  If  the  bubble  radius 
becomes  less  than  the  step  size  or  the  critical  radius,  the  bubble  vanishes. 

Derivation  of  Wash-out  Probability 

Our  intention  is  to  apply  the  model  to  study  gas  wash-out  from  animals  and  humans 
during  decompression.  Although  this  goal  requires  us  to  model  heterogeneously  distributed 
sources  and  sinks,  we  begin  with  tissue  taken  to  be  homogeneously  perfused  for  simplicity. 

This  assumption  implies  that  the  amount  of  gas  washed  out  per  unit  time  from  the 
tissue  is  inversely  proportional  to  the  tissue  time  constant,  T.  Thus  during  a  single  diffusion 
step  of  time  tstep,  tstep/T  of  the  material  will  wash  out.  This  allows  us  to  define  pwash  =  tstep/T 
as  the  probability  of  washing  out  on  a  single  step  so  that  we  can  calculate  the  wash-out 
probability  after  n^  steps  as, 
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Pcycle=PWashO-  +(1  -Pwash^1  ~PWash)2+-  +Pwash^  ^wastf^ 


(53) 


which  can  be  expressed  as, 


cycle  P  wash 


1-(1  -Pwash)^ 
1  “C1  -Pwash) 


(54) 


which  reduces  to, 


Pcycle^-Pwashf^  (55) 

Alternatively  for  ease  of  calculation,  this  can  be  expressed  as  an  exponential  since  if  e' 
pwash  is  expanded  as  a  Taylor  series  it  is  easy  to  see  that  for  small  values  of  pwash,  (l-pwash)  is 
approximately  equal  to  e'pwash.  In  fact,  if  pwash<  0.075  this  approximation  is  accurate  to  within 
4%.  Since  pwash  is  generally  smaller  than  this  we  can  rewrite  pcycle  as, 

p  =  1  -g  (56) 

cycle  c 

Since  ncycle-pwash  =  ncyclc-tstep/I  =  tcycle/T,  this  can  be  written  as, 


JjycU 


P  =l-e 

r cycle  1  e 


(57) 


The  inert  gas  concentration  in  tissue  after  wash-out  (Cadj)  can  be  obtained  by 
recognizing  that  Pcycle  is  simply  the  fraction  of  the  concentration  above  the  asymptotic  value 
(concentration  when  t  =  oo)  removed  during  wash-out,  so  that, 
or 
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p  =  CQ  Cadj 

cycle  C0-C„ 


Cadi  CQ  PcyCie(CQ  CJ 


(58) 

(59) 


where 

C0  =  concentration  before  adjusting  for  wash-out 
Co,  =  concentration  when  t=oo 


and  Co,  is  computed  from  the  arterial  inert  gas  tension,  (Painert),  which  we  assume  to  be  equal 
to  the  alveolar  inert  gas  tension,  (PAinert)  calculated  as  (22), 

<6#) 


where, 


%  -  <«> 

Combining  these  two  equations  yields, 

(«) 


where: 

PACo2  =  Alveolar  carbon  dioxide  partial  pressure  (mmHg) 
Pa02  =  Alveolar  oxygen  partial  pressure  (mmHg) 

Fi02  =  Inspired  oxygen  fraction 

so  that, 


p 
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(63) 
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and 
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Model  Testing 

We  tested  our  model  by  first  verifying  that  the  equilibrium  conditions  predicted  by 
Tikuisis  and  Ward  (17,18)  were  met  for  closed  systems  with  one  gas  as  calculated  by  the 
equations  in  Appendix  A.  Wash-out  was  then  simulated  without  the  presence  of  a  bubble  and 
compared  with  the  expected  single  exponential  curve.  Next,  we  compared  the  solution  with 
that  of  a  partial  differential  equation  model.  Finally,  we  developed  measures  for 
summarizing  the  inherently  variable  Monte  Carlo  simulations  and  examined  the  model's 
performance  at  different  levels  of  precision. 

For  comparison  of  the  time  course  of  bubble  evolution  with  that  produced  by  a 
standard  method,  we  solve  the  partial  differential  equation  (PDE): 

dPlk  =  +  2  (65) 

dt  dr2  r  dr 


where  Pliq  is  the  pressure  of  inert  gas  in  the  tissue,  D  is  the  diffusion  coefficient,  t  is  time, 
and  r  the  radial  coordinate.  The  boundary  at  the  bubble-liquid  interface  satisfies  the 
condition  that  the  inert  gas  partial  pressure  near  the  bubble  in  the  liquid  is  the  same  as  that  in 
the  bubble  and  also  satisfies  the  condition  that  the  flux  into  the  bubble  is  dependent  on  the 
gradient  in  the  liquid  near  the  bubble  so  that, 
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~dt 


(66) 


dP,. 

=  -aDA. — — 
b  dr 


' bub 


Where  n,  is  the  number  of  moles  in  the  tissue,  a  is  the  solubility,  and  Abub  is  the  surface  area 
of  the  bubble  of  radius  rbub.  The  outer  boundary  is  a  reflection  boundary.  The  internal 
bubble  pressure  is  given  by  Equation  3.  The  real  positive  root  of  Equation  52  is  the  new 
bubble  radius.  Gas  wash-out  is  assumed  to  be  uniform  throughout  the  liquid.  It  is  calculated 
assuming  the  rate  of  gas  elimination  is  proportional  to  the  amount  of  gas  present.  We  have 
developed  a  solution  to  this  system  using  the  Crank-Nicholson  method. 

Because  Monte  Carlo  models  are  based  on  random  samples  from  the  distributions  of 
interest,  no  two  simulations  will  be  the  same.  As  a  consequence  it  is  necessary  to  do 
multiple  runs  and  use  summary  measures  for  comparisons.  To  estimate  the  variability  in  our 
model  predictions  we  did  10  runs  at  3  different  levels  of  particle  density  precision:  2  x  108 
particles/cm3,  4  x  108  particles/cm3,  and  8  x  108  particles/cm3.  (For  comparison,  water  in 
equilibrium  with  nitrogen  gas  at  1  AT  A  and  room  temperature  has  approximately  5  x  1017 
molecules  of  N2/cm3.)  When  summarizing  the  10  runs  at  a  given  level  of  precision,  five 
measures  are  of  particular  interest:  the  maximum  radius  achieved  (rmM),  time  to  maximum 
radius  (t^),  time  to  bubble  dissolution  (tdis),  mean  transit  time  (TT),  and  relative  dispersion 
of  the  transit  times  (RD)  (2325).  RD  is  equal  to  the  standard  deviation  of  the  TT  divided  by 
the  mean  TT.  From  the  moment  of  the  step  change  in  P^,  t^  and  tdis  are  measured. 

Mean  TT  and  RD  do  not  represent  the  parameters  of  a  single  transit  time  distribution,  as  they 
do  under  normobaric  conditions.  This  is  because  the  evolving  bubble  continuously  changes 
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the  transit  time  distribution  of  the  system.  Nevertheless,  this  approach  was  selected  because 
these  measures  represent  a  simple  means  for  quantifying  the  effect  of  bubble  evolution  on  gas 
transit  through  tissue,  which  is  readily  comparable  with  prior  work. 

The  transit  time  probability  density  function  describes  the  distribution  of  the  transit 
times  of  all  the  inert  gas  particles  in  the  tissue.  If  we  express  the  transit  time  probability 
density  function  f(t)  as, 

M  <67) 

where  the  B;  and  $  constants,  the  wash-out  function  g(t)  is  given  by, 

t 

g(f )  =  QcJ,\-jfix)dx)  («*) 


which  reduces  to, 

*«)  -  <*„£.  P,P^'#8‘  <69> 

where  Qca  is  a  scaling  factor  attributable  to  the  blood  flow  and  wash-in  concentration. 
Furthermore,  we  know  that, 


TT  =  £.  P(P? 


(70) 
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and 


m  __  '/[2(E,  Wft-nl  (71) 

TT 

In  order  to  estimate  TT  and  RD,  we  fit  g(t)  to  the  simulated  wash-out  curves  produced  by  the 
model,  then  use  the  exponential  parameters  to  calculate  our  estimates  of  the  mean  transit  time 
and  relative  dispersion  of  transit  times  (26). 

Model  Parameters 

The  parameters  for  the  model  can  be  divided  into  three  categories:  (1)  the  physical 
variables,  (2)  the  boundary  values  immediately  after  the  bubble  is  formed,  and  (3)  the 
simulation  control  inputs.  The  physical  variables  include  the  pressure  profile,  surface 
tension,  bubble  density,  washout  time  constant,  the  diffusion  coefficient,  and  the  Ostwald 
solubility  coefficient  for  the  liquid.  The  immediate  post-bubble  formation  boundary 
conditions  can  range  from  uniform  gas  distribution  to  depletion  of  the  inner  most  shells  of  the 
inert  gas  in  the  bubble  but  maintaining  partial  pressure  equilibrium  across  the  boundary.  The 
simulation  control  parameters,  grid  size,  and  bubble  growth  cycle  time  were  selected  so  as  to 
assure  the  stability  of  the  solution.  The  simulation  results  were  considered  stable  when 
changes  in  either  variable  did  not  substantially  alter  the  bubble  evolution  curves. 

The  values  of  the  parameters  used  for  the  test  runs  are  summarized  in  Table  1  by 
equation  number  and  reference.  Those  values  marked  with  a  PDE,  were  used  in  the  Monte 
Carlo  -  PDE  comparison.  The  remainder  of  the  physiological  parameters  were  the  same  for 
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all  comparisons.  For  those  parameters  marked  "see  text"  in  the  reference  column,  the 
following  provides  the  rationale  for  the  values  selected. 

The  values  selected  for  Vn  were  chosen  for  computational  efficiency,  but  as  close  as 
possible  to  the  physiological  range.  Francis  et  al.  (31,32)  reported  bubble  densities  in  cross 
sections  of  dog  spinal  cords  ranging  from  0.007  ±  0.017  bubbles/mm2  to  0.251  ±  0.487 
bubbles/mm2.  This  corresponds  to  an  intrabubble  distance  of  approximately  10  mm  to  2  mm. 
If  these  values  are  assumed  to  also  represent  the  vertical  intrabubble  distance  and  bubbles  are 
assumed  to  be  uniformly  distributed,  the  bubble  density  would  range  from  1  bubble/cm3  to 
125  bubbles/cm3.  Since  some  cord  sections  had  1  bubble/mm2  the  maximum  density  may 
extend  to  1000  bubbles/cm3.  Low-density  systems  require  more  computer  execution  time 
because  of  the  larger  liquid  volume  they  represent.  Densities  in  the  physiological  range  are 
impractical  given  our  current  computer  technology.  Consequently,  we  arbitrarily  selected 
2000  bubbles/cm3  and  8000  bubbles/cm3  even  though  such  densities  might  exist  only  in 
extremely  severe  decompression  sickness.  These  represent  tissue  volumes  of  0.0005  and 
0.000125  cm3,  respectively. 

The  pressure  step  from  (2  ATA)  to  (1  ATA)  was  selected  arbitrarily. 

The  value  for  the  initial  bubble  size  (rbub)  for  the  equilibrium  test  was  6.5  microns. 

Ten  to  twenty  percent  of  the  bubbles  with  this  starting  radius  shrink  below  the  step  size 
because  of  random  fluctuations  before  they  are  able  to  reach  a  size  close  to  the  equilibrium 
radius.  We  increased  the  starting  radius  to  9.2  microns  for  the  Monte  Carlo  and  PDE 
comparison  tests  to  avoid  this  problem. 

Diffusion  coefficients  for  nitrogen  range  from  0.6  x  10'5  cm2/sec  to  2.2  x  10'5  cm2/sec 
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for  tissue  with  water  content  between  65%  and  100%.  Tendon  and  bone  are  approximately 
65%  water,  muscle  is  approximately  75%  water,  and  brain  is  approximately  85%  water  with 
corresponding  diffusion  coefficients  of  about  0.6  x  10‘5,  0.75  x  10'5,  and  1.5  x  10'5  cm2/sec 
(29,30).  The  choice  of  diffusion  coefficient  affects  computer  execution  speed,  i.e.,  in  that 
slower  diffusion  requires  less  computer  execution  time.  We  ran  the  equilibrium  test  at  0.5  x 
10'5  cm2/sec,  but  to  remain  in  a  physiological  range,  we  used  a  diffusion  coefficient  of  0.75  x 
10'5  cm2/sec  for  the  Monte  Carlo  -  PDE  comparison  tests. 

The  grid  size  for  all  tests  was  2.5  microns,  which  corresponds  to  a  diffusion  step  size, 
(a),  of  4.3  microns.  Diffusion  step  size  in  the  Monte  Carlo  model  and  the  grid  size  in  the 
PDE  model  were  selected  to  assure  that  smaller  values  would  not  change  the  solution. 

The  number  of  steps  in  a  bubble  growth  cycle  was  set  at  50  for  both  the  equilibrium 
and  Monte  Carlo  -  PDE  comparison  tests.  Because  the  diffusion  coefficients  were  different 
for  the  two  tests  while  the  grid  size  remained  the  same,  the  time  for  a  single  step  as 
calculated  by  Equation  26  was  different.  This  results  in  bubble  growth  cycle  time  of  0.3125 
seconds  for  the  equilibrium  test,  but  0.21  seconds  for  the  Monte  Carlo  -  PDE  comparison 
test. 

Longer  tissue  wash-out  times  require  more  computer  execution  time,  so  we  selected 
tissue  times  that  were  relatively  fast  but  close  to  the  physiologic  range.  This  resulted  in  our 
choice  of  2.5  and  10  minutes  for  the  tests. 

RESULTS 

The  critical  (rc)  and  equilibrium  (re)  radii  for  our  system  with  a  volume  of  0.000125 
cm3  (8,000  bubbles/cm3)  and  0.0005  cm3  (2,000  bubbles/cm3)  without  mass  transfer  were  1.4 


43 


microns  and  70.4  microns,  respectively.  Figure  9  demonstrates  the  model's  ability  to  reach 
the  predicted  equilibrium  radius  (re)  when  the  system  is  closed  to  mass  transfer.  A  bubble 
with  a  radius  greater  than  re  shrinks  to  re  and  a  bubble  smaller  than  re  grows  to  it. 

The  simulated  inert  gas  wash-out  without  a  bubble  precisely  matched  the  expected 
single  exponential  curve  predicted  by  a  well-stirred  model. 

The  Monte  Carlo  and  PDE  predictions  are  compared  in  Table  2  and  an  illustrative 
condition  is  presented  in  Figures  10  and  1 1 .  The  four  combinations  of  bubble  density  and 
wash-out  time  referred  to  as  conditions  I-IV  are  presented  in  Table  2.  Condition  IV  was  not 
calculated  for  the  PDE  model  because  of  the  computational  time  required.  Mean  transit  time 
and  relative  dispersion  for  condition  II  and  the  PDE  model  were  not  calculated  because  of 
difficulties  encountered  in  fitting  the  wash-out  curves  with  Equation  69.  Values  for  the 
Monte  Carlo  mean  rmax  compare  favorably  with  the  PDE  r^  for  conditions  I-m.  Because  of 
the  variability  in  the  Monte  Carlo  method,  the  remainder  of  the  measures  are  not  as  closely 
matched.  However,  almost  all  fall  within  one  standard  deviation  of  the  mean.  In  general, 
the  predictions  of  the  Monte  Carlo  and  PDE  models  are  similar  as  is  graphically  illustrated 
for  condition  I  in  Figures  10  and  11.  Figures  12  and  13  show  the  variability  of  the  model 
predictions  according  to  the  number  of  particles  in  each  simulation  for  bubble  radius  verses 
time;  Figure  12  has  2  x  108  particles/cm3  and  Figure  11  has  8  x  108  particles/cm3.  Figures  14 
and  15  are  the  corresponding  graphs  for  gas  wash-out.  They  show  the  number  of  moles  of 
gas  in  the  system  over  time.  Each  figure  contains  curves  from  10  separate  simulations  using 
the  same  starting  conditions.  The  expected  single  exponential  curves  for  inert  gas  wash-out 
without  the  presence  of  a  bubble  are  shown  for  comparison  in  the  bottom  of  Figures  14  and 
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FIGURE  9  -  Graph  of  simulation  output  for  bubble  evolution  in  a  system  closed  to  inert  gas 
transport  with  initial  bubble  radii  greater  than  and  less  than  the  expected  equilibrium  bubble 
radius.  Tissue  volume  is  1.25  x  10^  cm3.  Note  that  both  simulations  converge  to  the 
equilibrium  radius  predicted  by  the  quartic  polynomial  in  Appendix  A  (70.4  microns). 
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TABLE  2-  Comparison  of  Monte  Carlo  and  partial  differential  equation  model 
predictions  of  the  time  course  of  bubble  evolution 


Condition 

I 

II 

III 

IV 

Bubble  Density 
(bubbles/cm3) 

8,000 

8,000 

2,000 

2,000 

Wash-out  Time 

(sec) 

150 

600 

150 

600 

critical  radius 
(microns) 

1.4 

1.4 

1.4 

1.4 

equilibrium  radius 
(microns) 

70.4 

70.4 

112.3 

112.3 

Monte  Carlo  Model 

number  of  runs 

10 

10 

10 

10 

fmax1  (microns) 

43.0±2.7 

61.7±2.1 

42.8+5.3 

80.3  +  1.5 

tmax'  (Sec) 

342.1  ±146.1 

545.6±133.9 

282.7±  107.9 

820.8  +  143.6 

tdis1  (sec) 

950.4±  140.5 

3799.6±349.4 

771.3±203.3 

4918.5  +  845.8 

Mean  Transit  Time1,5 

286.8±37.23 

2039. 2  ± 227. 33 

177.1±17.43 

1403.0+ 164.74 

Relative  Dispersion1,5 

1.44±0.133 

1.16±0.043 

1.19±0.143 

1.0* 

Partial  Differential 
Equation  Model 

rmax  (microns) 

tmax  (sec) 

tdis  (sec) 

Mean  Transit 
Time23-5 

Relative 
Dispersion2,3,5 

1.  Mean  ±  s.d.,  coefficient  of  variation  (s.d./mean)).  2.  Estimated  from  1  (1,2,  and  1  runs  respectively)  and  2  exponential  fits  of 

simulated  wash-out  data.  3.  Mean  transit  time  for  liquid  without  a  bubble  is  150  seconds  with  a  relative  dispersion  of  1.  4.  Estimated  from 
1-exponential  fit  of  simulated  wash-out  data.  4.  Estimated  from  1-exponential  fit  of  simulated  wash-out  data.  5.  Mean  transit  time  for 
liquid  without  a  bubble  is  equal  to  the  wash-out  time.  The  relative  dispersion  is  1. Dissolution  for  the  Monte  Carlo  model  occurs  when  rbub 
is  4.3  microns.  *  -  not  studied. 


40.5 

59.5 

43.2 

* 

236.0 

504.3 

277.2 

* 

965.0 

3930.0 

1013.0 

* 

301.3+26.5 

* 

185.5+3.0 

* 

1.26+0.07 

* 

1.21+0.009 

* 
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FIGURE  13  -  Graph  of  bubble  radius  versus  time  for  higher  precision  simulation  (8  x  108 
particles/cm3)-  Tissue  volume  is  1.25  x  lO"4  cm3. 
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FIGURE  14  -  Graph  of  gas  wash-out  versus  time  for  low  precision  simulation  (2  x  108 
particles/cm3).  Thin  line  represents  the  expected  wash-out  for  the  tissue  without  a  bubble. 
Tissue  volume  is  1.25  x  104  cm3. 
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FIGURE  15  -  Graph  of  gas  wash-out  versus  time  for  higher  precision  simulation  (8  x  108 
particles/cm3).  Thin  line  represents  the  expected  wash-out  for  the  tissue  without  a  bubble. 
Tissue  volume  is  1.25x  10-4  cm3. 
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15.  These  figures  illustrate  the  variability  inherent  in  a  method  based  on  Monte  Carlo 
methods,  and  the  need  for  summary  measures  for  comparing  results. 

The  summary  measures,  mean  ±  S.D.  and  coefficients  of  variation  (S.D./mean)  of 
the  maximum  radius  achieved  (rmM),  time  to  maximum  radius  (tmax),  time  to  bubble 
dissolution  (tdis),  mean  transit  time,  and  relative  dispersion  of  the  transit  times,  for  10 
simulation  runs  at  each  of  the  three  particle  density  precisions  (2  x  108/cm3,  4  x  108/cm3,  and 
8  x  108/cm3)  are  presented  in  Table  3.  In  our  test  system,  10  -  20  %  of  the  bubbles  with  a 
starting  radius  of  6.5  microns  shrink  below  the  step  size  before  they  are  able  to  reach  a  size 
close  to  the  equilibrium  radius.  As  a  consequence,  the  results  are  biased  towards  fast¬ 
growing  bubbles,  although  the  size  of  the  bias  is  probably  small. 

These  results  illustrate  that  increasing  particle  density  in  a  bubble  growth  cycle  does 
not  affect  the  precision  of  outcome  measures  equally;  since  each  particle's  transit  through  the 
module  is  independent,  we  might  expect  that  the  precision  of  the  outcome  measures  improves 
proportionally  to  the  square  root  of  the  relative  increase  in  particle  density.  This  would 
imply  that  the  coefficient  of  variation  (cv)  for  the  highest  particle  density  run  in  Table  3 
would  be  half  the  cv  for  the  lowest  density  run.  While  r^,  t^  and  tdis  approximate  this 
expected  improvement,  the  mean  transit  time  and  relative  dispersion  of  the  transit  times  do 
not.  Since  precision  for  all  measures  will  improve  proportionally  with  the  square  root  of  the 
number  of  complete  runs,  multiple  runs  at  low  particle  density  may  be  a  more  efficient  use  of 
computer  resources  than  a  few  high  precision  runs  for  estimating  the  mean  transit  time  and 
the  relative  dispersion.  More  exhaustive  testing  will  be  required  to  determine  this 
definitively. 


53 


TABLE  3:  Precision  of  Outcome  Measures  as  a  Function  of  Particle  Density 


2  x  108/cm3 
(n=10) 

fmax1  (M) 

34.5 ±7.6,  0.22 

228.3 ±127.2,  0.56 

tdis^S) 

657.1  ±273.4,  041 

Mean  Transit  Time1,2,3 

204.3 ±45.3,  0.22 

Relative  Dispersion1 ,2,3 

1 .21  ±0.18,  0.15 

Particle  Density 

4  x  108/cm3 
(n=10) 

33.3±4.6,  0.14 

8  x  108/cm3 
(n=10) 

32.5±4.7,  0.14 

271.9±122.2,  0.45 

280.6±83.8,  0.30 

793.4±274.6,  0.35 

838.1±177.4,0.21 

196.9±25.4,  0.13 

200.4±28.9,  0.14 

1.23±0.17,  0.14 

1.27±0.17,  0.13 

1.  Mean  ±  s.d.,  coefficient  of  variation  (s.d./mean)). 

2.  Estimated  from  1  (1,2,  and  1  runs  respectively)  and  2  exponential  fits  of  simulated  wash-out  data. 

3.  Mean  transit  time  for  liquid  without  a  bubble  is  150  seconds  with  a  relative  dispersion  of  1.4.  Estimated  from 
1 -exponential  fit  of  simulated  wash-out  data. 

4.  Estimated  from  1-exponential  fit  of  simulated  wash-out  data. 

5.  Mean  transit  time  for  liquid  without  a  bubble  is  equal  to  the  wash-out  time.  The  relative  dispersion  is  1. 
Dissolution  for  the  Monte  Carlo  model  occurs  when  rbub  is  4.3  microns. 
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The  program  execution  speed  is  directly  proportional  to  the  number  of  particles 
simulated.  Since  larger  volumes  require  more  particles  to  obtain  the  same  number  of  bubble- 
liquid  transitions,  larger  volumes  require  longer  execution  speeds.  The  relationship  among 
the  number  of  steps  taken  in  a  bubble  growth  cycle,  the  grid  size  and  volume  and  execution 
speed  is  more  complex  since  all  these  factors  influence  the  proportion  of  particles  that  will  be 
able  to  take  a  bigstep  in  the  liquid  region. 

DISCUSSION 

Our  ultimate  goal  is  the  development  of  a  physiologically  based  method  of 
constructing  decompression  schedules.  Since  the  effect  of  extravascular  bubbles  on  inert  gas 
exchange  must  be  included  in  such  an  endeavor,  we  set  out  to  develop  a  method  for 
constructing  models  of  bubble  evolution  during  decompression  that  would  allow  for 
previously  excluded  details  of  the  physiological  environment  such  as  complex  microvascular 
architectures. 

As  a  first  step  in  this  direction,  this  report  introduces  a  model  of  bubble  evolution 
during  decompression  in  a  homogenous  liquid  based  on  a  Monte  Carlo  simulation  of 
diffusion.  The  major  advantage  of  this  method  is  that  in  theory,  details  of  the  biological 
environment  can  be  incorporated  that  are  not  feasible  using  traditional  modeling  methods. 

The  major  disadvantage  is  that  Monte  Carlo  simulations  require  large  amounts  of  computer 
execution  time.  Other  disadvantages  include  the  inherent  variability  in  the  model  predictions, 
which  requires  multiple  runs  and  the  use  of  summary  measures  for  comparisons,  and  the 
dependence  of  execution  time  on  the  volume  of  the  liquid  being  simulated.  Limitations  of  the 
current  simplified  model  are  the  assumptions  of  liquid  homogeneity,  uniform  gas  wash-out, 
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uniform  bubble  size  and  distribution,  radial  symmetry,  and  the  use  of  a  spherical  region 
which  is  not  space  filling. 

Two  features  of  this  approach  are  novel  to  Monte  Carlo  models  of  diffusion  and 
contribute  to  its  relative  efficiency.  First,  use  of  the  probability  distribution  of  the  distance  a 
particle  travels  in  the  liquid  phase  to  replace  many  single  diffusion  steps  dramatically 
increases  the  program's  execution  speed.  It  also  invites  the  development  of  other 
distributional  approaches.  Second,  treating  the  bubble  as  a  high  solubility  region  allows  us  to 
model  gas  particle  flux  without  developing  a  detailed  kinetic  model  of  particle  behavior  inside 
the  bubble. 

The  results  presented  in  this  report  confirm  the  fact  that  a  Monte  Carlo  approach  can 
produce  the  equilibrium  results  expected  from  independent  calculations  and  compares 
favorably  with  the  predicted  time  course  of  bubble  evolution  predicted  by  a  PDE  model. 

Future  model  development  could  progress  in  several  ways.  Effort  at  improving  the 
efficiency  of  the  algorithm  could  result  in  execution  times  that  are  competitive  with  traditional 
methods.  One  approach  is  to  attempt  to  generalize  the  probability  distribution  function  so  as 
to  calculate  in  advance  a  distribution  that  includes  both  diffusion  and  transition  into  and  out 
of  the  bubble.  Sampling  from  the  distribution  of  particle  residence  times  in  the  bubble  could 
also  improve  efficiency.  The  final  step  would  be  to  develop  a  single  distribution  which 
would  combine  all  of  these  elements.  For  the  special  case  of  the  spherically  symmetric 
bubble-liquid  module,  the  problem  can  be  reduced  to  a  one  dimensional  random  walk.  The 
theoretical  background  for  these  approaches  is  outlined  in  Appendix  D. 

Moving  the  program  to  a  computer  capable  of  parallel  or  vector  processing  might  also 
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improve  execution  speed  since  each  particle's  path  is  independent  of  every  other  particle  and 
could  be  calculated  separately  (33). 

In  the  current  implementation  we  treat  metabolic  gases  as  constants  for  simplicity. 
Alternatively,  each  type  of  gas  could  be  treated  stochastically  based  on  its  solubility,  tissue 
time  and  diffusion  coefficient  and  tracked  separately  through  a  bubble  growth  cycle.  An 
additional  "metabolic  conversion"  probability  would  need  to  be  applied  to  account  for  the 
consumption  of  oxygen  and  production  of  carbon  dioxide.  Multiple  inert  gases  could  also  be 
handled  individually.  Multiple  gas  simulations  are  restricted  by  execution  speed  since  each 
new  gas  requires  the  addition  of  an  equal  number  of  particles  to  maintain  precision.  So  a 
model  with  nitrogen,  helium  and  oxygen  would  require  three  times  as  many  particles  as  the 
simulations  in  this  report. 

The  influence  of  bubble  interactions  may  be  another  area  of  study.  This  would 
require  a  further  partitioning  of  the  liquid  region  in  order  to  more  precisely  record  the 
gradient  in  the  liquid.  Since  each  liquid  region  requires  a  minimum  number  of  particles  to 
ensure  a  large  enough  number  of  particle  transitions  between  regions  in  each  bubble  growth 
cycle,  this  approach  would  also  require  considerably  more  execution  time.  A  similar 
approach  could  be  taken  to  include  details  of  the  tissue  micro-architecture.  Incorporation  of  a 
symmetrical,  space-filling  liquid  region  would  allow  generalization  from  one  liquid  module  to 
an  entire  liquid  slab  without  concern  about  errors  introduced  by  the  use  of  non  space-filling 
liquid  regions  currently  employed  (34). 

Ultimate  model  testing  depends  on  experiments  conducted  in  in  vitro  and  in  vivo 
systems.  If  experimental  methods  with  sufficient  resolving  power  are  developed,  direct 
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testing  could  be  conducted  by  comparing  predicted  rmax,  tmax  and  tdis  with  experimental  results 
An  indirect  approach  could  be  taken  with  present  technology  by  comparing  predicted  transit 
time  measures  with  those  observed  in  animal  experiments  conducted  with  radiolabeled  gases 
(24,25). 
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APPENDIX  A:  Determination  of  Equilibrium  Radii 

The  polynomial  (17,18)  which  describes  the  2  equilibrium  radii  in  a  closed  system  with 
one  inert  gas  is: 


a4-r4  +  a3r3  +  axr  +  aQ  =  0  (Al) 


where: 


a4  =  4H-P(l-PJP)l3cs-k-T 
a3  =  SIL‘y/3‘Cs‘k‘T 
-  ni-p^p)-Nics 
a0  =  1TVIP 


and 

r  =  bubble  equilibrium  radii,  (cm) 

P  =  ambient  inert  gas  partial  pressure  after  pressure  step,  (mmHg) 

PVap  =  vapor  pressure  of  water,  (constant),  (mmHg) 

cs  =  saturation  concentration  of  inert  gas  in  the  solvent,  (moles/cm3) 

For  weak  solutions  cs  =  c,  P/KH,  where: 

c,  =  molar  concentration  of  the  solvent,  (moles/cm3) 

Kh  =  Henry's  coefficient  =  PH/([(k-T/(VL-a-PH))+ 1]'1),  (mmHg)  (35) 

where: 

VL  =  molar  volume  of  solvent  (1/c,),  cm3/mole 
a  =  Ostwald  solubility  coefficient 
P„  =  760  mmHg 

k  =  Universal  gas  constant  in  (cm3  mmHg)/(mole  K) 

T  =  temperature  in  Kelvin 
y  =  surface  tension,  (mmHg-cm) 

V  =  volume  of  depleted  region  (inverse  of  bubble  density),  (cm3) 

N  =  number  of  moles  of  inert  gas  in  the  volume  V  at  cs,  (moles) 
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APPENDIX  B:  Derivation  of  Initial  Particle  Placement  Coordinates 


We  need  to  create  the  desired  distribution  from  a  uniform  random  variable  between  0  and 
1  generated  by  the  computer.  Because  our  volume  is  a  sphere,  we  begin  in  spherical  coordinates 
p,  0  and  (|)  where  p  is  the  radial  coordinate,  0  the  azimuthal  angle  between  0  and  71,  and  (j) 
orthogonal  to  0,  between  0  and  271.  These  variables  are  mutually  independent,  therefore  we  must 
generate  a  distribution  for  each  of  them; 

Let  U  be  a  uniform  random  variable  on  [0, 1] ,  with  the  probability  density  function  (PDF) 
f(u)  and  the  cumulative  distribution  function  (CDF)  F(u),  where  u  is  a  particular  value  of  U. 
Let  p,  0  and  (j)  be  random  variables  with  PDF's  gp,  g0  and  g^  and  CDF's  Gp,  G0  and  G^  such 
that  p,  0  and  (j)  are  distributed  uniformly  by  volume.  For  illustrative  purposes  we  will 
demonstrate  the  procedure  in  detail  for  p,  although  similar  expressions  can  be  written  for  the 
other  variables,  0  and  (j).  By  definition, 


Profc{p0^p^p1}=Gp(p1)-Gp(p0)= 


volume  of  a  sphere  with  p^p^pj 
total  volume  of  the  sphere 


where  p,  and  p,  are  arbitrarily  selected  values. 

We  need  to  find  the  monotonic  transformations  hp(u),  h^u)  and  h^u),  such  that  p=hp(u), 
0=he(u)  and  (|)=h^(u).  So  for  example,  we  must  find  hp(u)  such  that, 


volume  ofaspherewith  p0^p^pj 
Prob{  P  o^W‘P.1-  roto,TOtoco/rtespAere 

^p(Pi-po)  (B2> 

3 
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where  R  is  the  radius  of  the  sphere. 

To  generate  the  desired  distribution  from  the  uniform  distribution,  we  require, 
Prob{p<hp(u)}  =  Prob{U<u}  which  by  definition  of  a  CDF  implies  that  Gp(hp(u))  =  F(u). 


Gp(*p(“»  =  F<“>  (B3a) 

Differentiating  both  sides  of  (B3a)  and  applying  the  chain  rule  gives, 

*p(p)A'p(“)l/W  (B3b) 

By  substituting  hp'l(p)  for  u  and  recognizing  that  f(u)=l  on  [0,1],  we  obtain, 


S.(P)=t; 


h'p(Ap‘‘(P» 


Integrating  gives, 


G„(P  )=/ - \ - dp 

A'plA,  (P» 

Combining  equations  Bl,  B2,  and  B5  for  pwe  obtain, 


Pi 


f - - - dp  = 

J  h\h-\ p) 


Po 


4  71  (  3  3X 
— (Pi"Po) 

471  Z>3 


Taking  the  derivative  with  respect  to  Pj  yields, 


1 


3  Pi 


hp(h-p\p )) 


(B4) 


(B5) 


(B6) 


(B7) 
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which  is  equivalent  to, 


1  _  3/tp(fQ  (B8) 

*>) "  *3 


since  hp(u)  =  pby  definition.  We  can  rearrange  to  obtain  the  differential  equation, 


dhp_  r*  (B9) 

du 


which  can  be  integrated  to  obtain  the  expression, 


1  (BIO) 

h=Ru  3 

p 


Recognizing  that  for  0, 


2n 


volume  of  a  sphere  with  O^O^Oj  3  ^  ^  cosfOp+cosfOo)) 


total  volume  of  the  sphere 


—R3 


(Bll) 


and  for  (J), 


volume  of  a  sphere  with  ^<jjx 
total  volume  of  the  sphere 


(B12) 


similar  approaches  can  be  taken  to  obtain  results  for  these  other  variables.  Solving  for  h  in  each 
case  we  obtain, 
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he(u)  =  0  =  arccos(l-2u) 


(B13) 


and 


h^u)  =  <t>  =  27 i-u 


(B14) 
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APPENDIX  C:  Approaches  to  Increasing  Execution  Speed 
1.  Distribution  of  Particle  Residence  Times  in  the  Bubble 

The  probability  of  the  particle  being  available  to  exit  the  bubble  after  n  steps  can  be 
calculated.  During  a  single  step  the  particle  has  one  chance  of  being  available  to  exit,  so  the 
probability  of  being  available  after  n  steps  can  be  calculated  as, 

„  ,  (Cl) 

Pn ^PeJ1  +V  ~ Pexit )  +d  ~Pexi)+- ^eJ-1  ~PexuT  ) 


which  can  be  expressed  as, 


pn=P 


exit 


l-d-Pex,/ 

I'd  -PeJ 


(C2) 


From  this  we  can  solve  for  (l-pexi,)n, 

Pexit 

then  applying  logarithms  we  obtain, 


(C3) 


Ind-Ad-d-p^))) 

Pexit _ 

^~Pexi) 


(C4) 


Since  pexit=l-(l-pexit), 


n_  ln(l-P„)  (C5) 

"  M1  ^exi) 

so  that  if  we  randomly  generate  Pn  between  0  and  1,  the  number  of  steps  before  the  particle  is 
available  can  be  calculated. 


70 


Alternatively,  for  pexit<D.075  the  geometric  distribution  can  be  approximated  by  the 
exponential  distribution  (21)  so  that, 

»=— —  -ln(w)  (C6) 

Petit 

where  u  is  a  uniform  random  variable. 

If  the  number  of  steps  generated  from  the  distribution  is  greater  than  the  number  of  steps 
in  a  bubble  growth  cycle,  n^,  the  particle  remains  unavailable  to  exit  during  the  bubble  growth 
cycle.  If  it  is  less  than  ncycle  then  for  (n^  -  n)  the  particle  takes  individual  steps. 

This  approach  can  be  applied  if  pexit  is  defined  as  described  in  the  main  body  of  the  text 
or  if  the  additional  pout  factor  were  included. 

2.  Distribution  of  particles  outside  the  bubble 

If  x  is  the  random  variable  that  describes  the  distance  outside  the  bubble  boundary,  we 
can  obtain  an  expression  for  the  distribution  of  the  particles  outside  the  bubble  by  integrating 
P(t  |  g)  with  respect  to  z  and  normalizing  the  result  by  the  expression  for  pexit, 


H[x)  = 


/ 


3  g  [2a(/?  +z)  -  (a 2  +2Rz  +z  2)]  [R  +z] 
4a/?3 

a  fl(12/?2-a2) 

16/?3 


dz 


(C7) 


which  reduces  to, 


jj,.  _  -x(3x2  Sax2  +  12Rx2 +l2xR2  -24axR+6a2x+12Ra2 -24aR2)  (C8) 

W=  a\\2R2-a2) 
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3.  One  dimensional  random  walk  for  a  spherically  symmetric  bubble-liquid  module 

The  algorithm  described  in  the  main  body  of  this  report  is  general  enough  to 
accommodate  a  heterogeneous,  asymmetric  diffusion  environment.  However,  a  spherically 
symmetric  bubble-liquid  module  can  be  studied  as  a  first  approximation  to  a  more  complicated 
geometry.  Because  of  the  symmetry  in  this  environment,  the  three-dimensional  random  walk  can 
be  reduced  to  a  one-dimensional  random  walk  in  order  to  increase  the  efficiency  of  the 
algorithm.  A  subroutine  FASTSTEP  has  been  developed  that  implements  this  one-dimensional 
algorithm. 

Subroutine  FASTSTEP  uses  the  law  of  cosines  to  generate  a  new  position  r'  for  a  particle 
initially  at  a  distance  r  from  the  origin,  and  stepping  a  distance  (a)  at  an  angle  (j)  to  the  radius 
vector,  as  shown  in  Figure  14.  The  equation  giving  r'  is  then, 

(C9) 

r'2  =  a2+r2-2’a7-cos(<|>) 

where  the  cosine  term  is  randomly  generated  from  a  uniform  random  distribution  by, 

cos(<j>)  =  1  -  2 -RAN {II)  (CIO) 

and  where  II  is  a  seed  for  the  random  number  generating  function  RAN(x). 

Using  the  geometry  of  Figure  15 ,  subroutine  FASTREFL  takes  particles  that  have  stepped 
out  of  the  system  and  return  them  back  into  the  system  using  reflection  from  the  outer  boundary. 
In  Figure  15,  a  particle  initially  at  point  P  steps  to  point  P',  which  is  at  a  distance  greater  than 
R,  the  radius  of  the  system.  In  order  to  reflect  the  particle  from  the  outer  boundary  to  point  P' ' , 
FASTREFL: 
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(1)  Uses  the  quantities  r,  R,  and  cos(<J))  (which  is  returned  from  subroutine 


(2) 


(3) 


FASTSTEP)  and  the  law  of  cosines  to  find  distance  dl5 

df  =  r2+R2-2-rR-cos(§) 

Uses  the  law  of  sines  to  determine  angle  0, 

sin(0)  =  —  sin(<J>) 

R 

Applies  the  law  of  cosines  a  final  time  to  determine  r", 
r"2  =  cg+R2-2-d2R-cos(Q) 


where  d2  =  a  -  dt. 


(Cll) 


(Cl  2) 


(Cl  3) 
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APPENDIX  D:  Fortran  Source  Code 


The  program  was  written  using  a  Fortran  compiler  for  Avalon  accelerator  boards  installed 
in  a  Vax  3800.  The  Fortran  source  code  for  the  program  is  listed  in  Section  Dl.  In  order  to 
execute  the  program,  in  addition  to  the  executable  code,  an  input  file  such  as  the  example  in 
Section  D2,  and  a  file  "DISTPROB.DAT"  (which  contains  the  values  of  the  diffusion  probability 
distribution  function  obtained  from  numerical  integration)  are  required.  This  file  is  produced  by 
the  program  listed  in  Section  D3.  An  additional  binary  file,  "BUBMOD.DMP",  must  also  exist 
prior  to  execution.  This  file  contains  information  on  the  gradient  that  can  be  used  to  restart  the 
program  in  the  middle  of  a  run.  The  statement,  iprofile=l  in  the  input  file  tells  the  program 
to  use  the  information  in  the  file.  The  output  file  created  by  the  program  contains  information 
on  the  time  course  of  the  bubble  radius.  A  program  listed  in  Section  D4  can  create  alternative 
initial  boundary  value  conditions  for  the  liquid  surrounding  the  bubble.  Section  D5  contains 
subroutines  that  can  be  used  to  implement  the  one-dimensional  random  walk  version  of  the 
program. 
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Section  Dl.  FORTRAN  source  cock  for  Monte  Carlo  bubble  simulation. 

C-  This  program  simulates  bubble  growth  beginning  with  a  bubble  of 
C-  critical  radius  as  defined  by  Tikuisis.  The  bubble  is  represented  as 
C-  a  region  of  high  solubility  surrounded  by  a  "depleted  region"  from 
C-  which  the  inert  gas  is  extracted.  These  two  regions  are  placed  in 
C-  a  third  region  which  represents  the  surrounding  tissue. 


implicit  real  *8  (a-h,o-z) 

common  /coml/xmin,xmax,ymin,ymax,zmin,zmax,rbub,rdep 
common/com2/pi,il 
common/distarr/distance(10000) 
common/sincos/fact,  scale ,  sinmat(0 : 9999)  ,cosmat(0 : 9999) 
common/com3/IA(3 ,24)  ,aa(3 ,24) 
common  /com5/twopi,onethird 
real*8  N2,N2dep,k,kH,molezone(500) 

Dimension  xx(3) ,  rzonebeg(0 : 500)  ,rzoneend(0 : 500) ,  widezone(500) 
dimension  strtzone(500),endzone(500),ffaczone(500),conczone(500) 
dimension  cummsum(0:500),vzone(500),rzsq(0:500) 
data  IA/1,1, 1,1, 1,-1, 1,-1, 1,1, -1,-1, -1,1, 1,-1, 1,-1, -1,-1, 

+  1,-1, -1,-1, 

+  1,1, 1,-1, 1,1, 1,-1, 1,-1, -1,1, 1,1, -1,-1, 1,-1, 1,-1, -1,-1, -1,-1, 

+  1,1,1, 1,1, -1,-1, 1,1, -1,1, -1,1, -1,1,1, -1,-1, -1,-1, 1,-1, -1,-1/ 

C-  Define  constants 


zed=0.d0 

Pi =4.d0*datan(l  .dO) 
twopi=2.d0*pi 
scale =1 0000.  dO 
fact=scale/twopi 
call  gensine 


c 

c  change  12/16/92  to  include  oxygen  and  co2  as  additional  gas 
c  components  of  the  bubble,  using  values  from  van  liew,  pg.  336 
c 

c  ph2o  =  6.2  kPa  =  46.52  mmHg 

c  po2  =5.3  kPa  =  39.76  mmHg 

c  pco2  =  5.9  kPa  =  44.26  mmHg 

c  pco21ung  =  5.3  kPa  =  39.76  mmHg 
c  pvap=  46.52d0  +  39.76d0  +  44.26d0 

c  read  in  from  input  file  below  1/18/93 
c  ph2o=46.52d0 

c  po2=39.76d0 

c  pco2=44.26d0 

c  pco21ung = 39 . 76d0 
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c  pvap=ph2o+po2+pco2 
C-  vapor  pressure  of  water  in  bubble  @  37  C  in  mmHg 
C1125  =0.055346 
C11=0. 055140 

C-  pure  solvent  cone  of  water  @  37  C  in  moles/cc 
k=62358.0 

C-  Boltzman's  gas  constant  in  mmHG-cc/mole-K 
C-  based  on  PV=kNT  and  STP  conditions 
R=62358.0 

C-  Universal  gas  constant  in  mmHg-cc/mole-K 
C-  based  on  published  value  of  0.08205  L-atm/mole-K 
T=37. 0+273.1 
C-  temperature  in  Kelvin 
onethird=l.d0/3.d0 
C  Pvap=Pvap25 

C  cll=cll25 

C-  Input  initial  conditions  and  control  variables. 

C-  ngrow= number  of  bubble  growth  cycles, 

C-  ntrip= number  of  independent  trips  thru  tissue,  nstep= number  of  steps  per 
C-  trip,  Surten= surface  tension  of  bubble  (dynes/cm),  PCbub= partition 
C-  coefficient  of  bubble  (1/ostwald  coefficient),  D=diffusion  coefficient 
C-  (cm  sq/sec),  Voldep0= initial  volume  of  the  depleted  region  (cu  cc), 

C-  Voltiss  =  volume  of  cubical  tissue  region 

C-  P01= initial  tissue  pressure  (mmHg),  Pstep= decompression  step  (mmHg) 
Rewind(l) 

Read( 1 , 900)ngrow ,  ntrip ,  nstep ,  surten,  Pcbub ,  D ,  gr  idsize ,  VoldepO , 

+  Voltiss,P01,Pstep,bubfact 

900  Format(I12/I12/I12/gl5.7/gl5.7/gl5.7/gl5.7/gl5.7/ 

+  g15.7/gl5.7/gl5.7/gl5.7) 

read(l ,  1900)iprofile,nmicro,washoutt,pctn2 

1900  format(il2/il2/gl5.7/gl5.7) 

read(l ,  1901)switch,volref,fracsd,ph2o,po2,pco2,pco21ung 

1901  format(gl5.7) 
close(l) 

pvap =ph2o +po2 +pco2 

C-  Input  the  array  of  diffusion  distances  as  a  function  of  probability 
C-  given  in  terms  of  the  standard  deviation 

Open(unit=  1  ,file  =  'distprob.daf  ,type  =  'old') 

Read(l  ,899)(distance(i),i  =  1 , 10000) 

899  Format(fl7.2) 

Close(l) 
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vol = voltiss + voldepO 
C-  Calculate  dimensions  of  cubical  region 
call  resize(vol,cubeside) 

C-  adjust  units  to  mmHG-cm  knowing  760  mmHg  =  1.013X10^6  dyn/cm^2 
Surten = Surten/ 1332.8 

C-  Calculate  Henry's  coefficient  from  PCbub 

KH = 760 .  /((((R*T)/ (( 1  ./cll)*(l  ./PCbub)  *760.))  + 1)**-1 .) 

C-  Calculate  stepsize 
c  modify  stepping  array 

call  regrid(ia,aa,gridsize , stepsize) 

C-  Calculate  bubble  growth  cycle  time  and  tissue  washout  time 
Dtime = stepsize*stepsize/(6 .  *D) 
eyelet = nstep*Dtime 
poofprob = 1  .dO-dexp(-cyclet/washoutt) 
if(  washoutt .  gt . 86400 .  d0)poofprob = 0 .  dO 
C-  Calculate  the  maximum  distance  a  particle  can  travel  by  diffusion 
onesd = dsqrt(2 .  d0*d*cyclet) 
diffdist = 5 .  dO  *onesd 
c 

c  figure  out  the  width  of  the  diffusion  shells,  switch  =  0  is  for 
c  a  dimensionless  run,  with  reference  volume  volref.  switch  <  >  0  is 
c  an  absolute  run,  with  parameters  unsealed. 

c 

if(switch.eq.0.d0)then 
width = onesd  *fracsd*(vol/volref)**onethird 
if(width.  gt.  onesd)width = onesd 
else 

width = onesd*ff  acsd 
end  if 

C-  Write  out  the  starting  conditions 

Open(unit=2,file  =  'for002.dat'  ,type  =  'new') 

Write(2,901)  Ngrow,Ntrip,Nstep,Surten*1332.8,PCbub,D, 

+  gridsize ,  stepsize ,  P01 ,  Pstep ,  V  oldepO ,  voltis  s 
Close(2) 

901  Format(3I9/ 5G 1 5 .4/4G  15.4) 

C-  Begin  a  bubble  growth  cycle 
time=0. 

P1=P01 

p  1  n2 = (p01-ph2o-pco21ung)  *pctn2 

pfh2 = (p01-pstep-ph2o-pco21ung)*pctn2 

Voldep= VoldepO 

il  =  182361  +2*secnds(0.0) 

lstrt=l 
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if(iprofile .  ne .  0)then 

c  for  a  nonuniform  initial  concentration  profile,  put  code  in  here 

open(unit=7,file=  'bubmod.dmp'  ,type  =  'old'  ,form=  'unformatted') 
rewind(7) 

read(7)lstrt,il,nzonebeg,time,volbub,calcpc,rbub,rdep,rbub0 

read(7)re,rc,currmole 

read(7)(fraczone(i),molezone(i),conczone(i),i=  1  ,nzonebeg) 

read(7)(rzonebeg(i),rzsq(i),vzone(i),i = 1  ,nzonebeg) 

close(7) 

lstrt = lstrt  + 1 

cs=cll*pln2/kh 

N2=Cs*Vol 

N2dep = Cs*Voldep 

orgmoles=N2 

If  (VoldepO.gt.vol)  orgmoles=n2dep 
Pl=Pl-Pstep 
csequil=cll*pfh2/kh 
rzonebeg(0)=0.d0 
rzoneend(O) = 0 .  dO 
bubmoles =molezone(l) 
open(unit=3,file  =  'for003.dat', type  =  'new') 
write(3,3000)0,time,rbub,volbub,calcpc, 

+  (fraczone(i),rzonebeg(i),conczone(i),i = 1  ,nzonebeg) 

close(3) 

depmoles = currmole-bubmoles 
bubfrac = bubmoles/currmole 
depff  ac = depmoles/currmole 
end  if 

Do  10  L= lstrt, ngrow 

C-  Calculate  number  of  moles  of  inert  gas  (N2)  at  start  of  the  cycle 
C-  in  the  entire  tissue  and  in  the  depleted  region  alone  (N2dep). 

If((L.eq. l).and.(iprofile.eq.O))  then 
5432  continue 

cs=cll*pln2/kh 
N2=Cs*Vol 
N2dep = Cs*Voldep 
orgmoles=N2 

If  (VoldepO.gt.vol)  orgmoles=n2dep 
currmole = orgmoles 
Pl=Pl-Pstep 
csequil =cll*pfh2/kh 
C-  Calculate  bubble  dimensions 
C-  First,  calculate  the  critical  and  equilibrium  radii 
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aO = 2.  *surten*Voldep/Pl 
Cs=Cll*Pl/KH 

al  =  Voldep*(l  .-Pvap/Pl)-N2dep/Cs 

a3 = 8 .  *pi*surten/  (3 .  *Cs  *k*T) 

a4 = 4 .  *pi*Pl  *(  1 .  -Pvap/Pl)/(3 .  *Cs  *k*T) 

c0=a0/a4 

cl=al/a4 

c3=a3/a4 

Call  roots(cO,cl,c3,rc,re) 
if  (rc.lt.O.)  goto  999 
C-  Assign  radius  to  test  bubble 

rbub = re*bubfact + rc  *(  1 .  dO-bubfact) 
rbub=dmaxl(rbub,(l  .5d0*stepsize)) 
rbubO=rbub 

VolBub =4.  *pi*rbub*rbub*rbub/3 . 
rdep = (3 .  *(Voldep + Volbub)/(4*pi))**(l  .12,.) 
c  check  here  to  make  sure  the  bubble  fits  inside  the  cube,  if  not,  spread 
c  the  excess  volume  onto  the  depleted  region,  and  make  it  larger 
if((rdep .  gt .  xmax) .  and .  (voltiss .  gt .  zed))then 
voldepO=vol 
voldep=voldepO 
voltiss = zed 
goto  5432 
end  if 

rzonebeg(0)=0.d0 
rzoneend(0) = 0 .  dO 

call  genzone(rbub,rdep,width,nmicro,rzonebeg,rzsq,vzone,nzonebeg) 
Bubmoles = Volbub*(Pl-Pvap +2.  *surten/rbub)/(R*T) 

Depmoles = N2Dep-Bubmoles 
Tismoles = N2-N2dep 
Bubfrac = Bubmoles/orgmoles 
Depff ac = Depmoles/orgmoles 
Tisfrac = Tismoles/orgmoles 
if(iprofile.eq.O)then 
fraczone(l) =bubffac 
molezone(  1 ) = bubmoles 
conczone(l) =bubmoles/volbub 
do  1100  izone=2,nzonebeg 

ffaczone(izone) = depffac*vzone(izone)/voldep 
molezone(izone) = fraczone(izone)*orgmoles 
conczone(izone) = depmoles/voldep 
1100  continue 
end  if 
End  if 
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call  genzone(rbub,rdep, width, nmicro,rzoneend,rzsq,vzone,nzoneend) 

C-  Adjust  the  dimensions  of  the  module 
Volcube = vol + volbub 
call  resize(volcube,cubeside) 

C-  Begin  the  random  walk  thru  the  module 
C-  Calculate  the  probability  of  staying  in  the  bubble 
shellsize = rbub-stepsize 
stepsiz2 = stepsize*stepsize 
shelsiz3 = shellsize*shellsize*shellsize 
rbub2 = rbub*rbub 
rbub3 = rbub2*rbub 
rdep2=rdep*rdep 

Pf=  1  ,dO/Pcbub*(rbub3-shelsiz3)/rbub3 
If  (rbub.lt.stepsize)  pf=l.dO/PCbub 
xldpf=-l.dO/pf 
Startbub=0. 

Startdep=0. 

Starttis=0. 

EndInbub=0. 

EndIndep=0. 

EndIntis=0. 

do  1005  izone=l,nzoneend 
strtzone(izone) =0.d0 
endzone(izone) = 0 .  dO 
1005  continue 

cummsum(O) = 0 .  dO 
do  1004  izone  =  l,nzonebeg 

cummsum(izone) = cummsum(izone- 1) + fraczone(izone) 

1004  continue 

tisffac=  1  ,dO-cummsum(nzonebeg) 
if(tisfrac.lt.0.dO)tisfrac =0.d0 
tismoles =tisffac*currmole 
C-  Calculate  the  number  of  particles  to  be  place 
Do  20  H=l,Ntrip 

C-  Place  a  particle  randomly  in  the  module  with  probability  weighted 
C-  by  mole  fraction. 

prob=ran(il) 

c 

c  new  code  to  place  particles  in  the  various  zones 
if(prob .  ge .  cummsum(nzonebeg))then 

call  putinbox(xx(l),xx(2),xx(3),zed,zed,zed,cubeside,partdist) 
else 
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do  1001  izone  =  l,nzonebeg 
if(prob .  It .  cummsum(izone))then 

call  zone(xx(l),xx(2),xx(3),zed,zed,zed,rzonebeg(izone-l), 
+  rzonebeg(izone),partdist) 

goto  1002 
end  if 

1001  continue 

1002  continue 
end  if 

C-  If  the  particle  is  placed  so  far  from  the  bubble  that  it  cannot 

C-  reach  it  in  nstep  steps  and  the  module  only  consists  of  the 

C-  bubble  and  a  depleted  region,  the  particle  will  stay  in  the  depleted 

C-  region. 

x=xx(l) 

y=xx(2) 

z=xx(3) 

If  (partdist.  gt.  (diffdist + rbub))  then 
Call  Bigstep(x,y,z,D,cyclet) 

C-  Check  to  see  if  the  particle  is  outside  the  spherical  depleted  region. 

C-  If  it  is  reflect  it  back. 

if(voltiss .  le .  zed)then 
else 

call  inbox(x,y,z) 
end  if 
Istep= nstep 

C-  The  particle  can  reach  the  bubble  by  diffusion 
Else  If  (partdist.gt.rbub)  then 
Istep = int((partdist-rbub)/stepsize) 

Steptime = Istep*Dtime 
Call  Bigstep(x,y,z,D, steptime) 
if(voltiss.le.zed)then 
else 

call  inbox(x,y,z) 
end  if 

C-  The  particle  is  in  the  bubble 
Else 
Istep =0 
End  if 

C-  Carry  out  random  walk  for  NSTEP  steps  of  "stepsize" 

Do  30  JJ = Istep +1, nstep 
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Inbub=0 

C-  Check  to  see  if  particle  is  in  the  bubble  by  comparing  the 
C-  particle  distance  from  origin  with  bubbble  radius 
if(distfun2(x,y,z,zed,zed,zed).le.rbub2) 

+  inbub=l 

C-  If  the  particle  is  in  the  bubble  place  in  the  transition  shell 
C-  based  on  the  probability  factor. 

If  (Inbub.eq.l)  then 
if  (ran(il).lt.pf)then 

if(rbub .  le .  stepsize)then 
call  putinorb(x ,y,z,zed,zed, zed , rbub) 
else 

call  zone(x,y,z, zed, zed, zed, rbub-stepsize, rbub, dummyr) 
end  if 

call  step(x,y,z,gridsize) 
end  if 
else 

C-  Take  a  step 

call  step(x,y,z,gridsize) 
c 

c  assume  now  that  if  particle  has  just  left  the  bubble,  it  cannot 
c  reach  the  boundaries  of  the  system  during  a  single  time  cycle, 
c  we  can  thus  eliminate  the  rest  of  this  do  loop. 

c 

End  if 

30  Continue 

If  (distfun2(x,y,z,zed,zed,zed).gt.rdep2) 

+  Call  Reflect(x,y,z,xx(l),xx(2),xx(3),zed,zed,zed,rdep) 

C-  Reassign  coordinates 
xx(l)=x 
xx(2)=y 
xx(3)=z 


C-  Keep  track  of  the  region  in  which  the  particle  ends 
dist = distfun2(x,  y  ,z ,  zed ,  zed ,  zed) 
do  1010  izone=l,nzoneend 
if(dist.  It.  rzsq(izone))then 

endzone(izone) =endzone(izone)  + 1  .d0 
goto  1011 
end  if 


82 


1010  continue 

1011  continue 

If  (dist.le.rbub2)  then 
Endlnbub = Endlnbub + 1 . 
else  If(dist.le.rdep2)  then 
Endlndep = Endlndep  + 1 . 
else 

Endlntis = Endlntis  + 1 . 
end  if 

20  Continue 

if((voltiss .  le .  zed) .  and .  (endintis .  gt .  zed))then 
open(unit =2,  file  =  'for002.dat'  ,type  =  'old’ , access  =  'append') 
write(2 , 8888)endintis ,  1 

8888  formate  accounting  error... ',fl0.2,' particles  in  tissue'/ 

+  '  on  the  ',i8,'th  cycle') 

stop 
end  if 

C-  Calculate  the  gas  content  in  each  region 
If  (L.eq.  1)  then 

open(unit =2,file  =  'for002.dat'  ,type = 'old' , access = 'append') 
write(2 ,904)time  ,rbub ,  volbub ,  voldep , 

+  (vol-voldep) 

write(2,904)cyclet,poofprob,washoutt 
Write  (2,905)  rc,re,rbub0 
close(2) 
end  if 
c 

c  this  is  where  washout  takes  place,  in  older  versions  (commented  out) 
c  washout  was  done  on  a  zone  by  zone  basis,  this  version  (2/23/93)  does 
c  washout  over  the  entire  depletion  region  at  the  same  time,  this  is  the 
c  same  mathematically  as  doing  them  separately, 
c 

dfnt = dfloat(ntrip) 
dfnt0=dfnt 
actual =endindep 
if(poofprob .  gt .  0 .  d0)then 
const = dfnt*csequil/currmole 
expected = const*voldep 

actual=endindep-(endindep-expected)*poofprob 
dfnt = actual + endinbub 
1025  continue 
1015  continue 
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if(voltiss .  gt .  zed)dfnt = dfnt + endintis 
currmole = currmole*dfht/dfntO 
end  if 

c 

c  end  of  washout  code 

c 

c  bookkeeping,  note  if  statement  to  get  proper  amount  in  bubble 
c 

do  1020  izone  =  l,nzoneend 

fraczone(izone) = endzone(izone)*actual/endindep/dfnt 
if(izone.eq.l)fraczone(izone)=endzone(izone)/dfnt 
molezone(izone) = ffaczone(izone)*currmole 
conczone(izone) = molezone(izone)/vzone(izone) 

1020  continue 

calcpc = ff  aczone(  1 )  *vzone(2)/(fraczone(2)  * vzone(  1 )) 
Bubfrac = EndlnBub/dfnt 
depfrac = actual/dfnt 
Tisffac = Endintis/dfnt 
bubmoles =currmole*bubfrac 
depmoles = currmole*depfrac 
tismoles = currmole*tisfrac 

C-  If  Bubfrac  is  very  small  it  will  soon  vanish  so  trap  it 
If  (bubfrac. lt.0. 0000000001)  goto  999 

C-  Calculate  the  new  volume  of  each  region 
C-  the  volume  of  the  depleted  region  plus  the  volume  of 
C-  the  supersaturated  tissue  region  is  a  constant,  so  bubble 
C-  growth  only  adds  to  max  and  min  dimensions,  but  not  volumes 
A3  =  (4.  *pi*(Pl-Pvap))/3 . 

A2 = (8 .  *pi*surten)/3 . 
aO = -molezone(l)*R*T 
CO = AO/ A3 
C2=A2/A3 

Call  Roots2(c0,c2,rbub) 
if(rbub.lt.stepsize)goto  999 
VolBub = (4*pi*(rbub)**3.)/3 . 
tis  VNew = T  ismoles/(C  1 1*P01/KHX 
Voldep = Vol-tisvnew 

If  (VoldepO.ge.vol)  Voldep =voldepO-tisvnew 
rdep = (3 .  *(Voldep + Volbub)/(4*pi))**(l  ./3 .) 
call  expand(rbub ,  rdep  .rzoneend ,  rzonebeg  ,nzoneend) 
nzonebeg = nzoneend 
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V  olcube = vol + volbub 
call  resize(volcube,cubeside) 

C-  Output  information  for  1  bubble  growth  cycle 
time = time  +  eyelet 

open(unit=2,file  =  'for002.dat'  ,type  =  'old' , access  =  'append') 

Wr ite(2 , 904)  time ,  rbub ,  ealepe ,  currmole 
close(2) 

904  Format(5G15.7) 

3000  formate  cycle  #  \i5,4gl2.4/50(3gl5.6/)) 

open(unit=7,file=  'bubmod.dmp' , type  =  'old'  ,form=  'unformatted') 
rewind(7) 

wr  ite(7)L ,  i  1 ,  nzonebeg ,  time ,  volbub ,  ealepe ,  rbub ,  rdep ,  rbubO 
write(7)re  ,rc , currmole 

write(7)(ffaczone(i),molezone(i),conczone(i),i =  1 , nzonebeg) 
write(7)(rzonebeg(i),rzsq(i),vzone(i),i = 1  .nzonebeg) 
close(unit=7) 

10  Continue 
999  Continue 

905  Format(3G15.4) 


Stop 

End 

c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

c  generate  sin  and  cos  matrices 
c 

subroutine  gensine 
implicit  real*8  (a-h,o-z) 

common/sincos/fact,  scale ,  sinmat(0 : 9999)  ,cosmat(0 : 9999) 
do  10  i =0,9999 
cosmat(i) = dcos(dfloat(i)/fact) 
sinmat(i) = dsin(dfloat(i)/fact) 

10  continue 
return 
end 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

subroutine  expand(rbub,rdep,rzl,rz2,nz) 
implicit  real  *8  (a-h,o-z) 


85 


dimension  rzl(0:500),rz2(0:500) 
common  /com2/pi,il 
common  /com5/twopi,onethird 
data  nzmax,eps/50, 1  .d-05/ 
coeff = 4 .  dO*pi  *onethird 
rz2(l)=rbub 
do  10  i=2,nz 

rz2(i)=(rzl(i)**3-rzl(i-l)**3+rz2(i-l)**3)**onethird 
10  continue 

do  20  i=nz+l,nzmax 
rz2(i)=0.d0 
20  continue 

if(dabs(rz2(nz)-rdep) .  gt.  eps)then 
write(6, 100)nz,rz2(nz),rdep 

100  format('  expansion  error:  rz2(',i2,')  =  ',fl0.3, 

+  '  rdep  =  \fl0.3) 

stop 
end  if 
return 
end 
c 
c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

subroutine  genzone(rbub ,  rdep ,  width ,  nm ,  rzone ,  rzsq , vzone , nzone) 
implicit  real*8  (a-h,o-z) 
dimension  rzone(0 : 500)  ,rzsq(0 : 500) ,  vzone(500) 
common/com2/pi,  i  1 
common  /com5/twopi,onethird 
nzone  =  1 
rzone(0)=0.d0 
rzone(l)=rbub 
coeff = 4 .  d0*pi*onethird 
vzone(l) =coeff*rzone(l)**3 
vtot=vzone(l) 
do  10  i= 2,500 
nzone = nzone +1 
rzone(i) = rzone(i- 1 ) + width 
if(rzone(i) .  gt.rdep)rzone(i) = rdep 
v = coeff*rzone(i)**3 
vzone(i)=v-vtot 
if(rzone(i).ge.rdep)goto  20 
vtot=v 

10  continue 
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20  continue 

do  40  i=0,nzone 

rzsq(i) = rzone(i)*rzone(i) 

40  continue 

return 
end 
c 

£******************************************************************** 
Subroutine  Roots(c0,cl,c3,rc,re) 
implicit  real*8  (a-h,o-z) 

C-  This  subroutine  calculates  the  roots  of  the  Tikuisis  polynomial. 

C-  First,  find  the  minimum  because  1  root  is  below  it  and  the  other 
C-  root  is  above  it. 

100  format('  #  \i5) 

If  (cl.ge.0.)  then 
rc=-999. 
else 

xl=-c0/cl 
10  x=xl 

g=cl+x*x*(3.*c3+4.*x) 
dgdx = x*(6 .  *c3  + 12 .  *x) 
xl=x-g/dgdx 

if  ((abs(g).gt.0.0000001).or.(abs(xl-x).gt.0.0000001)) 

+  goto  10 

C-  Now  apply  Newton's  method  to  find  the  roots 
f =rcfunc(xl  ,c0,cl  ,c3) 
if  (f.gt.0.)  then 
rc=-998 
else 

do  1  j  =  l,2 
if  (j  eq.l)  then 
mult=-l 
else 

mult=l 
end  if 

x2=xl*l.+0.1*mult 
20  x=x2 

fO = rcfunc(x,c0,c  1  ,c3) 
fl=cl+x*x*(3.*c3+4.*x) 
dxl=-f0/fl 
x2=dxl+x 

if  (abs(l.-(x2/x)).gt.0.0001)  goto  20 
if  (mult.lt.0)  rc=x2 
if  (mult.gt.0)  re=x2 
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1  continue 

end  if 
end  if 
return 
end 


Real*8  Function  rcfimc(x,cO, cl, c3) 

implicit  real*8  (a-h,o-z) 

z = cO + x*(c  1  + x*x*(c3 + x)) 

rcfunc=z 

return 

end 


***  *  *  *  *  *  *  **********  *  *********  *********  *  *  *  **  ****  *  ****  ****  ***  *  ****  *  ****  *  * 


Subroutine  Roots2(c0,c2,rbub) 
implicit  real*8  (a-h,o-z) 

C-  This  subroutine  calculates  the  roots  of  the  polynomial  describing 
C-  bubble  radius. 

C-  Apply  Newton's  method  to  find  the  roots 
x2=rbub 


20  x=x2 

fO = rbubfunc(x,c0,c2) 
fl=x*(2.*c2+3.*x) 
dxl=-fO/fl 
x2=dxl+x 

if  (abs(l.-(x2/x)).gt.0.0001)  goto  20 
rbub=x2 
return 
end 


Real*8  Function  rbubfunc(x,c0,c2) 
implicit  real*8  (a-h,o-z) 
z = cO + x*x*(c2 + x) 
rbubfunc=z 
return 
end 
c 
c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  subroutine  to  reflect  a  particle  from  the  inside  surface  of  a  sphere 
c 

c  (x,y,z)  current  position  of  the  particle 

c  (xo,yo,zo)  previous  position  of  particle 
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c  (xc,yc,zc)  center  of  the  bubble/depleted  region 

c  r  radius  of  the  depleted  region 

c 

subroutine  reflect(x,y,z,xo,yo,zo,xc,yc,zc,r) 
implicit  real*8  (a-h,o-z) 
iperm=0 
xosave=xo 
yosave=yo 
zosave=zo 
10  continue 

delx=x-xo 
dely=y-yo 
delz=z-zo 

if((delx .  eq .  0 .  dO) .  and .  (dely .  eq .  0 .  dO) .  and .  (delz .  eq .  0 .  d0))retum 
if(delz.eq.0.d0)then 
call  permute(x,y,z,xo,yo,zo,xc,yc,zc) 
iperm=iperm+l 
goto  10 
end  if 

cx=delx/delz 
cy= dely  /delz 
dx=xo-xc-cx*zo 
dy=yo-yc-cy*zo 
aq = cx*cx +cy*cy + 1  .dO 
bq = 2 .  d0*(cx*dx + cy  *dy-zc) 
cq = dx*dx + dy  *dy + zc  *zc-r  *r 
root = dsqrt(bq*bq-4 .  d0*aq*cq) 
zplus = (-bq + root)/(2 .  d0*aq) 
zminus = (-bq-root)/(2.d0*aq) 
xplus = cx*(zplus-zo) + xo 
yplus = cy  *(zplus-zo) + yo 
xminus = cx*(zminus-zo) + xo 
yminus = cy  *(zminus-zo) + yo 
dplus = distfun2(x,y  ,z , xplus , yplus ,  zplus) 
dminus = distfun2(x,y  ,z, xminus , yminus, zminus) 
if(dplus .  It .  dminus)then 
xp= xplus 
yp= yplus 
zp= zplus 
else 

xp= xminus 
yp= yminus 
zp= zminus 
end  if 
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cpx=xp-xc 

cpy=yp-yc 

cpz=zp-zc 

ptx=x-xp 

pty=y-yp 

ptz=z-zp 

c  get  dot  product  between  these  two  vectors 
dot = cpx*ptx + cpy  *pty + cpz  *ptz 
dt2or2 =2.d0*dot/(r*r) 
x=x-dt2or2*cpx 
y=y-dt2or2*cpy 
z=z-dt2or2*cpz 

if(distfunc(x,y,z,xc,yc,zc).gt.r)then 
xo=xp 
yo=yp 
zo=zp 
goto  10 
end  if 

if(iperm.ne.O)then 

do  20  i=jmod(iperm+l,3),jniod(iperm+3,3) 
call  permute(x,y,z,xo,yo,zo,xc,yc,zc) 

20  continue 
end  if 
xo=xosave 
yo=yosave 
zo=zosave 
return 
end 
c 
c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

subroutine  permute(al ,bl ,cl ,a2,b2,c2,a3 ,b3,c3) 

implicit  real*8  (a-h,o-z) 

savel=al 

save2=a2 

save3=a3 

al=bl 

bl=cl 

cl=savel 

a2=b2 

b2=c2 

c2=save2 

a3=b3 
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b3=c3 

c3=save3 

return 

end 

c 

c 

c 

c  subroutine  to  take  a  random  step  to  a  solid  angle  on  a  sphere 
c  of  radius  a 

c 

subroutine  Stepout(x,y,z,x0,y0,z0,bigr2,a) 
implicit  real*8  (a-h,o-z) 
common  /com2/pi,il 

common/sincos/fact,  scale ,  sinmat(0 : 9999) ,  cosmat(0 : 9999) 

c 

c  get  the  distance  from  the  center  of  the  sphere 


rx=x-xO 

ry=y-yO 

rz=z-zO 

d2 = rx*rx + ry  *ry + rz*rz 
d=dsqrt(d2) 

costhetb = (bigr2-a*a-d2)/(2  ,dO*a*d) 

ranmax=(l  .d0-costhetb)/2.d0 

costheta = (1  .d0-2.d0*ranmax*ran(il)) 

sintheta =dsqrt(l  .dO-costheta*costheta) 

delr=a*costheta 

delperp = a*sintheta 

urx=rx/d 

ury=ry/d 

urz=rz/d 

rperp =dsqrt(urx*urx +ury*ury) 
if(rperp.gt.0.d0)then 
ulx=-ury/rperp 
uly=urx/rperp 
else 

ulx=l.dO 
uly=0.d0 
end  if 

u2x=uly*urz 
u2y=-ulx*urz 
u2z =ulx*ury-urx*uly 
index = scale  *ran(i  1 ) 


91 


anorm=delr/d 

x = x + delperp*(cosmat(index)*ulx + sinmat(index)*u2x) + anorm*rx 
y = y +delperp*(cosmat(index)*uly + sinmat(index)*u2y) + anorm*ry 
z = z + delperp*(sinmat(index)*u2z) + anorm*rz 
return 
end 
c 
c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

c  subroutine  to  move  distance  a  from  a  point  in  the  boundary  region 
c  away  from  the  center  of  a  sphere, 
c 

subroutine  outstep(x,y,z,xO,yO,zO,a) 

implicit  real*8  (a-h,o-z) 

dx=x-xO 

dy=y-yO 

dz=z-zO 

anorm = a/dsqrt(dx*dx  4-  dy  *dy + dz*dz) 
x=x+anorm*dx 
y=y+anorm*dy 
z=z+anorm*dz 
return 
end 
c 
c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

c  subroutine  to  take  a  single  large  diffusion  step  in  three 
c  dimensions,  uses  subroutine  stepit,  the  random  direction 

c  stepping  routine,  and  also  calculates  a  normally  distributed 

c  distance  that  the  particle  steps, 
c 

subroutine  bigstep(x,y,z,diffcoef,t) 
implicit  real*8  (a-h,o-z) 
common  /com2/pi,il 
common/distarr/distance  ( 1 0000) 
c 

c  the  distance  stepped  as  a  function  of  the  probability,  given  in  terms 
c  of  the  standard  deviation,  is  located  infile  distprob.dat 
c 

index=aint(l  .d0+ 10000. dO*ran(il)) 
sd =dsqrt(2.d0*diffcoef*t) 
d = sd*distance(index) 
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call  stepit(x,y,z,d) 

return 

end 

c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

c  subroutine  to  take  a  random  step  of  length  rO  in  3  dimensions 
c 

subroutine  stepit(x,y,z,rO) 
implicit  real*8  (a-h,o-z) 
common  /com2/pi,il 

common/sincos/fact,  scale ,  sinmat(0 : 9999) ,  cosmat(0: 9999) 
common  /com5/twopi,onethird 
xr  =  1  .d0-2.d0*ran(il) 
sintheta = dsqrt(  1 .  dO-xr  *xr) 
index = scale  *ran(i  1 ) 
x = x + rO*sintheta*cosmat(index) 
y = y + rO*sintheta*sinmat(index) 
z=z+rO*xr 
return 
end 
c 
c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

c  subroutine  to  position  a  particle  randomly  within  a  spherical 
c  shell  with  outer  radius  router, and  inner  radius  rinner 

c 

c 

subroutine  zone(x,y,z,xO,yO,zO,rinner,router,rho) 
implicit  real  *8  (a-h,o-z) 
common  /com2/pi,il 

common/sincos/fact ,  scale ,  sinmat(0 : 9999)  ,cosmat(0 : 9999) 
common  /com5/twopi,onethird 
c 

c  use  spherical  coordinates  rho, theta, phi 
c 

x=xO 

y=yO 

z=zO 

if  (rinner .  le .  0 .  dO)then 
rho=0.d0 
return 
end  if 
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rin3 =rinner*rinner*rinner 
rout3 = router  *router*router 
range = rout3-rin3 

rho = (rin3  +  range  *ran(i  1 ))  *  *onethird 

call  stepit(x,y,z,rho) 

return 

end 


c 

c  subroutine  to  position  a  particle  randomly,  within  a  spherical 
c  shell  with  outerradius  bigr,and  inner  radius  bigr-a 
c 
c 

subroutine  shell(x,y,z,xO,yO,zO,rO,a) 
implicit  real*8  (a-h,o-z) 
common  /com2/pi,il 

common/sincos/fact,  scale ,  sinmat(0 : 9999)  ,cosmat(0 : 9999) 
common  /com5/twopi,onethird 

c 

c  use  spherical  coordinates  rho,theta,phi 

c 

if(r0.1t.a)then 

call  putinorb(x,y,z,xO,yO,zO,rO) 
return 
end  if 
bma=rO-a 

bma3 =bma*bma*bma 
b3=r0*r0*r0 
range =b3-bma3 

rho = (bma3 + range  *ran(i  1 ))  *  *onethird 

phi = twopi*ran(il) 

theta =dacos(l  .d0-2.d0*ran(il)) 

x = xO + rho*dsin(theta)  *dcos(phi) 

y = yO + rho*dsin(theta)*dsin(phi) 

z = zO + rho  *dcos(theta) 

return 

end 


c 

C 


subroutine  step(x,y,z,grid) 
implicit  real  *8  (a-h,o-z) 
common/com2/p  i ,  i  1 
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common/com3/IA(3 ,24),aa(3 ,24) 

J  A  =  1 + int(24 .  *ran(i  1 )) 
x=x+aa(l,ja) 
y=y+aa(2,ja) 
z=z+aa(3,ja) 
return 
end 
c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

c  subroutine  to  test  for  presence  of  a  particle  in  a  sphere  of  radius  r 
c  centered  at  (xO,yO,zO),  given  position  (x,y,z) 
c 

subroutine  insphere(x,y,z,xO,yO,zO,r,in) 
implicit  real*8  (a-h,o-z) 
logical  in 
r2=r*r 

if(distfun2(x,y,z,x0,y0,z0).le.r2)then 

in=.true. 

else 

in  =.  false, 
end  if 
return 
end 

c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

c  subroutine  to  place  particle  randomly  within  a  sphere  of  radius  r 
c  centered  at  (xO,yO,zO) 

c 

subroutine  putinorb(x,y,z;xO,yO,zO,rO) 
implicit  real*8  (a-h,o-z) 
common  /com2/pi,il 

common/sincos/fact,  scale ,  sinmat(0 : 9999)  ,cosmat(0 : 9999) 
common  /com5/twopi,onethird 
c 

c  use  spherical  coordinates  rho, theta, phi 

c 

rho = rO*(ran(il)**onethird) 

x=xO 

y=yO 

z=zO 
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call  stepit(x,y,z,rho) 

return 

end 

c 

c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

c  subroutine  to  see  if  particle  has  left  the  box  and  if  so,  to 
c  reflect  it  from  the  boundary, 
c 

subroutine  inbox(x,y,z) 
implicit  real*8  (a-h,o-z) 

common  /coml/xmin,xmax,ymin,ymax,zmin,zmax,rbub,rdep 
c 

c  reflection  conditions 

c 

if(x .  It .  xmin)x = xmin + (xmin-x) 
if(x .  gt .  xmax)x = xmax + (xmax-x) 
if(y  .lt.ymin)y =ymin+ (ymin-y) 
if(y.gt.ymax)y =ymax+(ymax-y) 
if(z .  It .  zmin)z = zmin + (zmin-z) 
if(z .  gt .  zmax)z = zmax + (zmax-z) 
c 

c  translation  conditions 

c 

c  if(xx(l).lt.xmin)xx(l) =xmax-(xmin-xx(l)) 

c  if(xx(l).gt.xmax)xx(l) =xmin-(xmax-xx(l)) 

c  if(xx(2).lt.ymin)xx(2) =ymax-(ymin-xx(2)) 

c  if(xx(2) .  gt .  ymax)xx(2) = ymin-(ymin-xx(2)) 

c  if(xx(3) .  It .  zmin)xx(3) = zmax-(zmin-xx(3)) 

c  if(xx(3).gt.zmax)xx(3) =zmin-(zmax-xx(3)) 

c 

return 

end 


c 

c 


^st:*******************************!^*********  ******************************** 


c  function  to  find  the  distance  between  two  points 
c 


function  distfunc(x,y,z,xO,yO,zO) 
implicit  real*8  (a-h,o-z) 
dx=x-xO 


dy=y-yO 

dz=z-zO 
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distfunc = dsqrt(dx  *dx + dy  *dy + dz  *dz) 

return 

end 

c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

c  function  to  calculate  the  square  of  the  distance  between  points 

c 

function  distfun2(x,y,z,x0,y0,z0) 

implicit  real*8  (a-h,o-z) 

dx=x-xO 

dy=y-yO 

dz=z-zO 

distfun2 = dx*dx + dy  *dy + dz*dz 

return 

end 

c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

subroutine  regrid(ia , aa,gridsize , stepsize) 
implicit  real*8  (a-h,o-z) 
dimension  ia(3,24),aa(3,24) 
do  1901  iii=l,3 
do  1901  iiii = 1 ,24 

aa(iii,  iiii) = dfloat(ia(iii ,  iiii))  *gridsize 
1901  continue 

stepsize =gridsize*dsqrt(3  .d0) 
return 
end 
c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

subroutine  resize(vol,s) 
implicit  real*8  (a-h,o-z) 

common  /coml/xmin,xmax,ymin,ymax,zmin,zmax,rbub,rdep 
common  /com5/twopi,onethird 

s = (vol**onethird) 
xmax=s/2. 
xmin=-xmax 
ymax=xmax 
ymin=xmin 
zmax=xmax 
zmin=xmin 
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return 

end 

c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

subroutine  putinbox(x,y,z,xO,yO,zO,s,dist) 
implicit  real*8  (a-h,o-z) 

common  /coml/xmin,xmax,ymin,ymax,zmin,zmax,rbub,rdep 
common  /com2/pi,il 
10  continue 

x = xmin + s  *r  an(i  1 ) 
y = y  min + s  *r  an(i  1 ) 
z = zmin + s  *ran(i  1 ) 
dist=distfunc(x,y,z,xO,yO,zO) 
if(dist.lt.rdep)goto  10 
return 
end 
c 
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Section  D2.  Sample  input  file  for  source  code  of  Section  Dl. 


8000, 

number  of  cycles 

400000, 

number  of  trips  per  cycle 

50, 

number  of  steps  per  trip 

50.00, 

surface  tension 

70.00, 

partition  coefficient 

.0000050, 

diffusion  coeff 

.00025,  grid  size 

.  125E-03 ,  depletion  region  volume 

0.0,  extra-depletion  volume 

1520.00,  initial  inert  gas  partial  pressure 

760.00,  pressure  step 

.1706,  initial  bubble  radius  (as  fraction  of  equil.  radius) 

0,  iprofile  #  of  initial  concentration  profile 

0,  nmicro  if  of  sub-zones  to  be  used 


150., 

washoutt 

tissue  time  constant  in  seconds 

0.7900, 

pctn2 

percent  nitrogen  during  washout 

l.dO, 

switch  =0 

=  >  dimensionless  run 

.100e-02, 

volref 

reference  vol  for  dimless  run 

l.dO, 

fracsd 

fraction  of  onesd  for  shell  widths 

46.52d0, 

ph2o 

partial  pressure  water  vapor 

39.76d0, 

po2 

partial  pressure  oxygen  in  bubble 

44.26d0, 

pco2 

partial  pressure  C02  in  bubble 

39.76d0, 

pco21ung 

part,  press.  C02  in  lungs 

l.dO, 

bndfact 

boundary  factor  for  PDE3 

1/19/93  test  run  for  .ch9 
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Section  D3.  FORTRAN  source  code  for  generation  of  cunumdative  distribution  function  for  diffusion  in  three 
dimensions. 


c  program  to  generate  the  file  'distprob.dat',  which  contains  the 
c  number  of  standard  deviations  stepped  by  a  diffusing  particle  as 

c  a  function  of  probability,  it  uses  integration  routine  from  the 

c  numerical  recipes  directory  on  piggy, 
implicit  real*8  (a-h,o-z) 
npts  =  10000 
nstep= 1000000 
dfnp=dfloat(npts) 
sum=0.d0 
sumnext=0.d0 
xlow=0.d0 
delx =0.0001 

open(unit=  1, file  =  'distprob.dat'  ,type  =  'new') 
do  10  i=l,npts-l 
prob = dfloat(i)/dfnp 
xhigh=xlow 
do  20  j  =  l,nstep 
xhigh = xhigh + delx 
call  qrombfxlow, xhigh, sumnext) 
if((sum+sumnext).lt.prob)goto  20 
sum = sum + sumnext 
write(l ,  100)(xlow +xhigh)/2.d0,i 
100  format(fl5.6,il0) 

xlow= xhigh 
goto  10 

20  continue 
10  continue 
close(l) 
stop 
end 
c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

function  func(t) 
implicit  real*8  (a-h,o-z) 
data  init/0/ 
if(init.eq.0)then 
pi = 4 .  dO  *datan(  1 .  dO) 
const = dsqrt(2 .  dO/pi) 
init=l 


100 


end  if 

func = const*t*t*exp(-t*t/2 .  dO) 

return 

end 

c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

SUBROUTINE  QROMB(A,B,SS) 
implicit  real*8  (a-h,o-z) 

PARAMETER(EPS = 1  .E-6,  JMAX = 20JMAXP = JM  AX  + 1  ,K = 5  ,KM = 4) 
DIMENSION  S(JM AXP)  ,H(JM AXP) 

H(l)  =  l. 

DO  11  J=1,JMAX 
CALL  TRAPZD(A,B,S(J),J) 

IF  (J.GE.K)  THEN 
L=J-KM 

CALL  POLINT(H(L),S(L),K,0.,SS,DSS) 

IF  (ABS(DSS).LT.EPS*ABS(SS))  RETURN 
ENDIF 
S(J+1)=S(J) 

H(J + 1)=0.25*H(J) 

11  CONTINUE 

PAUSE  'Too  many  steps.' 

END 

c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

SUBROUTINE  TRAPZD(A,B,S,N) 
implicit  real  *8  (a-h,o-z) 

IF  (N.EQ.l)  THEN 
S = 0. 5  *(B-A)*(FUNC(A) + FUNC(B)) 

IT=1 
ELSE 
TNM=IT 
DEL = (B- A)/TNM 
X=A+0.5*DEL 
SUM=0. 

DO  11  J=1,IT 
SUM=SUM+FUNC(X) 

X=X+DEL 
11  CONTINUE 

S = 0 . 5  *(S + (B- A)  *SUM/TNM) 

IT =2*IT 
ENDIF 
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RETURN 

END 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

SUBROUTINE  POLINT(XA,YA,N,X,Y,DY) 
implicit  real*8  (a-h,o-z) 

PARAMETER  (NMAX=10) 

DIMENSION  X A(N) ,  Y A(N) ,  C(NM AX) ,  D(NM AX) 

NS  =  1 

DIF= ABS(X-XA(1)> 

DO  11  1=1, N 
DIFT = ABS(X-XA(I)) 

IF  (DIFT. LT. DIF)  THEN 
NS=I 
DIF = DIFT 
ENDIF 
C(I)=YA(I) 

D(I)=YA(I) 

11  CONTINUE 
Y=YA(NS) 

NS=NS-1 

DO  13  M  =  1,N-1 
DO  12  I=1,N-M 
HO=XA(I)-X 
HP = XA(I + M)-X 

w=ca+D-Da) 

DEN=HO-HP 

IF(DEN.EQ.O.)PAUSE 

DEN=W/DEN 

D(I)=HP*DEN 

C(I)=HO*DEN 

12  CONTINUE 

IF  (2*NS.LT.N-M)THEN 
DY=C(NS  +  1) 

ELSE 

DY=D(NS) 

NS=NS-1 

ENDIF 

Y=Y+DY 

13  CONTINUE 
RETURN 
END 
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Section  D4.  FORTRAN  source  code  for  the  generation  of  non-uniform  initial  gas  concentration  distributions. 


c  this  program  is  a  modification  of  the  beginning  of  the  bubble  simulation 
c  program,  it  is  designed  to  generate  a  particular  non-uniform  concentration 
c  gradient  for  use  by  the  main  program,  the  main  program  assumes  an  initial 
c  uniform  concentration  profile  wherein  the  gas  used  to  create  the  bubble 
c  is  drawn  from  throughout  the  tissue  volume,  this  program  generates  a 
c  profile  from  the  opposite  extreme;  the  gas  for  the  bubble  is  taken  from 
c  the  innermost  shells,  leading  to  a  step  function  profile,  shells  near  the 
c  bubble  have  initial  concentration  cs  =  concentration  of  gas  in  the  bubble, 
c  shells  far  from  the  bubble  have  concentration  csO,  in  equilibrium  with 
c  the  initial  pressure  of  inert  gas.  there  is  one  transition  shell  with  a 
c  concentration  between  these  two  values,  the  concentration  profile  is  then 
c  written  out  to  file  bubmod.dmp,  where  it  can  be  used  by  the  main  program 
c  by  setting  iprofile  =  1 
c 


implicit  real*8  (a-h,o-z) 

common  /coml/xmin,xmax,ymin,ymax,zmin,zmax,rbub,rdep 
common/com2/pi,il 
common/com3/IA(3 ,24),aa(3 ,24) 
common  /com5/twopi,onethird 
real*8  N2,N2dep,k,kH,molezone(500) 

Dimension  xx(3),rzonebeg(0:500),rzoneend(0:500),widezone(500) 
dimension  strtzone(500),endzone(500),ffaczone(500),conczone(500) 
dimension  cummsum(0 : 500) , vzone(500) ,  rzsq(0 : 500) 
data  IA/1,1, 1,1, 1,-1, 1,-1, 1,1, -1,-1, -1,1, 1,-1, 1,-1, -1,-1, 

+  1,-1, -1,-1, 

+  1,1, 1,-1,1,1,1,-1,1, -1,-1, 1,1, 1,-1, -1,1, -1,1, -1,-1,-1,-1,-1, 

+  1,1, 1,1, 1,-1, -1,1, 1,-1, 1,-1, 1,-1, 1,1, -1,-1, -1,-1, 1,-1, -1,-1/ 

C-  Define  constants 


zed=0.d0 

Pi=4.d0*datan(l  .dO) 
twopi=2.d0*pi 
fourpio3  =  4 .  d0*pi/3 .  dO 
onethird  =  1  .dO/3  .dO 
c  Pvap25 =23.756 

c  Pvap=47.067 

c  change  12/16/92  to  include  oxygen  and  co2  as  additional  gas 
c  components  of  the  bubble,  using  values  from  van  liew,  pg.  336 
c 

c  ph2o  =  6.2  kPa  =  46.52  mmHg 

c  po2  =  5.3  kPa  =  39.76  mmHg 

c  pco2  =5.9  kPa  =  44.26  mmHg 
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c  pco21ung  =  5.3  kPa  =  39.76  mmHg 
c  pvap=  46.52d0  +  39.76dO  +  44.26d0 
ph2o=46.52d0 
po2=39.76d0 
pco2  =44.26d0 
pco21ung = 39 . 76d0 
pvap = ph2o + po2 + pco2 

C-  vapor  pressure  of  water  in  bubble  @  37  C  in  mmHg 
C1125 =0.055346 
01=0.055140 

C-  pure  solvent  cone  of  water  @  37  C  in  moles/cc 
k=62358.0 

C-  Boltzman's  gas  constant  in  mmHG-cc/mole-K 
C-  based  on  PV=kNT  and  STP  conditions 
R=62358.0 

C-  Universal  gas  constant  in  mmHg-cc/mole-K 
C-  based  on  published  value  of  0.08205  L-atm/mole-K 
T=37.0+273.1 
C-  temperature  in  Kelvin 
C  Pvap=Pvap25 

C  cll=cll25 

C-  Input  initial  conditions  and  control  variables. 

C-  ngrow= number  of  bubble  growth  cycles, 

C-  ntrip= number  of  independent  trips  thru  tissue,  nstep= number  of  steps  per 
C-  trip,  Surten= surface  tension  of  bubble  (dynes/cm),  PCbub= partition 
C-  coefficient  of  bubble  (1/ostwald  coefficient),  D= diffusion  coefficient 
C-  (cm  sq/sec),  VoldepO= initial  volume  of  the  depleted  region  (cu  cc), 

C-  Voltiss  =  volume  of  cubical  tissue  region 

C-  P01= initial  tissue  pressure  (mmHg),  Pstep= decompression  step  (mmHg) 
Rewind(l) 

Read(l,900)ngrow,ntrip,nstep,surten,Pcbub,D,gridsize,Voldep0, 

+  Voltiss, P01,Pstep,bubfact 

900  Format(I9/I9/I9/F9.2/F9.2/F9.7/F9.7/F9.7/F9.7/F9.2/F9.2/F9.7) 

read(l ,  1900)iprofile,nmicro,washoutt,pctn2 
1900  format(il2/il2/gl5.7/gl5.7) 
close(l) 

vol = voltiss + voldepO 
C-  Calculate  dimensions  of  cubical  region 

C-  adjust  units  to  mmHG-cm  knowing  760  mmHg  =  1.013X10*6  dyn/cmA2 
Surten = Surten/ 1332.8 

C-  Calculate  Henry's  coefficient  from  PCbub 

KH =760./((((R*T)/((l  ./cll)*(l  ./PCbub)  *760.)) +!)**-!.) 
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C-  Calculate  stepsize 
c  modify  stepping  array 

call  regrid(ia,aa,gridsize, stepsize) 

C-  Calculate  bubble  growth  cycle  time  and  tissue  washout  time 
Dtime = stepsize*stepsize/(6 .  *D) 
cy  clet = nstep  *Dtime 

C-  Calculate  the  maximum  distance  a  particle  can  travel  by  diffusion 
onesd = dsqrt(2 .  dO  *d*cyclet) 

C-  Begin  a  bubble  growth  cycle 
time=0. 

P1=P01 

p  1  n2 = (p01-ph2o-pco21ung)  *pctn2 
pfn2 = (p01-pstep-ph2o-pco21ung)*pctn2 
Voldep=VoldepO 
lstrt  =  1 

C-  Calculate  number  of  moles  of  inert  gas  (N2)  at  start  of  the  cycle 
C-  in  the  entire  tissue  and  in  the  depleted  region  alone  (N2dep). 

5432  continue 

c  Cs=Cll*Pl/KH 

cs=cll*pln2/kh 
csO=cs 
N2=Cs*Vol 
N2dep = Cs*Voldep 
orgmoles=N2 

If  (VoldepO.gt.vol)  orgmoles=n2dep 
currmole = orgmoles 
Pl=Pl-Pstep 
csequil=cll*pfn2/kh 
C-  Calculate  bubble  dimensions 
C-  First,  calculate  the  critical  and  equilibrium  radii 
aO = 2 .  *surten*  Voldep/Pl 
c  Cs=Cll*Pl/KH 

cs = c  ll*(pl-pvap)/kh 
al  = Voldep*(l.-Pvap/Pl)-N2dep/Cs 
a3 = 8.  *pi*surten/(3 .  *Cs*k*T) 
a4 = 4 .  *pi*Pl*(  1 .  -Pvap/Pl)/(3 .  *Cs  *k*T) 
c0=a0/a4 
cl=al/a4 
c3=a3/a4 

Call  roots(c0,cl,c3,rc,re) 
if  (rc.lt.O.)  goto  999 
C-  Assign  radius  to  test  bubble 
rbub=re*bubfact 
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rbub = dmax  1  (rbub ,  (1 . 5dO*stepsize)) 
rbubO=rbub 

VolBub =4.  *pi*rbub*rbub*rbub/3 . 
c  rdep = (3 .  *(V  oldep + Volbub)/(4*pi))  *  *(1 . /3 . ) 

rdep = ((voldep + volbub)/fourpio3)  *  *onethird 
c  check  here  to  make  sure  the  bubble  fits  inside  the  cube,  if  not,  spread 
c  the  excess  volume  onto  the  depleted  region,  and  make  it  larger 
if((rdep .  gt .  xmax) .  and .  (voltiss .  gt .  zed))then 
voldepO=vol 
voldep =voldepO 
voltiss = zed 
goto  5432 
end  if 

rzonebeg(0)=0.d0 
rzoneend(0) = 0 .  dO 

call  genzone(rbub,rdep,onesd,nmicro,rzonebeg,rzsq,vzone,nzonebeg) 
Bubmoles = Volbub*(Pl-Pvap + 2 .  *surten/rbub)/(R*T) 

Depmoles = N2Dep-Bubmoles 
Tismoles = N2-N2dep 
Bubffac =Bubmoles/orgmoles 
Depfr  ac = Depmoles/orgmoles 
Tisffac = Tismoles/orgmoles 
ffaczone(l) =bubffac 
molezone(l) =bubmoles 
conczone(l) =bubmoles/volbub 
deles =csO-cs 
vinner = bubmoles/delcs 
vrunold=0.d0 
do  1100  izone=2,nzonebeg 
vrun = vrunold + vzone(izone) 
if(vrun.  le .  vinner)then 
conczone(izone) = cs 

else  if((vrunold .  le .  vinner) .  and .  (vinner .  It .  vrun))then 

conczone(izone) = ((vinner- vrunold)  *cs + (vrun-vinner)*csO)/ 

+  vzone(izone) 

else 

conczone(izone) = csO 
end  if 

vrunold = vrun 

molezone(izone) = conczone(izone)*vzone(izone) 
ff aczone(izone) = molezone(izone)/orgmoles 
1100  continue 

3000  formate  cycle  #  ’,i5,4gl2.4/50(3gl5.6/)) 

1=0 
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rbubO=rbub 

time=0.d0 

il  =237645 +2*secnds(0.0) 

calcpc = ff  aczone( 1 )  *vzone(2)/(fraczone(2)  * vzone(  1 )) 
open(unit=7,file  =  'bubmod.dmp'  ,type  =  'new' , form  =  'unformatted') 
rewind(7) 

write(7)L,il,nzonebeg, time, volbub, calcpc, rbub,rdep,rbubO 
write(7)re,rc,currmole 

write(7)(ffaczone(i),molezone(i),conczone(i),i  =  1  ,nzonebeg) 
write(7)(rzonebeg(i),rzsq(i),vzone(i),i  =  1  ,nzonebeg) 
close(unit=7) 

999  Continue 

Stop 

End 


c 

c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

subroutine  genzone(rbub,rdep,width,nm,rzone,rzsq,vzone,nzone) 
implicit  real*8  (a-h,o-z) 
dimension  rzone(0:500),rzsq(0:500),vzone(500) 
common/com2/pi ,  i  1 
common  /com5/twopi,onethird 
nzone = 1 
rzone(0)=0.d0 
rzone(l)=rbub 
coeff = 4  .dO*pi*onethird 
vzone(  1 ) = coeff*rzone(  1 )  *  *3 
vtot=vzone(l) 
do  10  i=2,50 
nzone = nzone +1 
rzone(i) = rzone(i- 1 ) + width 
if(rzone(i).gt.rdep)rzone(i)=rdep 
v =coeff*rzone(i)**3 
vzone(i)=v-vtot 
if(rzone(i).ge.rdep)goto  20 
vtot=v 

10  continue 

20  continue 

do  40  i=0, nzone 
rzsq(i) =rzone(i)*rzone(i) 
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40 


continue 

return 

end 


c 


Subroutine  Roots(c0,cl,c3,rc,re) 
implicit  real*8  (a-h,o-z) 

C-  This  subroutine  calculates  the  roots  of  the  Tikuisis  polynomial. 
C-  First,  find  the  minimum  because  1  root  is  below  it  and  the  other 
C-  root  is  above  it. 

100  format('  #  '  ,i5) 

If  (cl.ge.O.)  then 
rc=-999. 


else 

xl=-c0/cl 


10  x=xl 

g = c  1  +  x*x*(3 .  *c3  +4 .  *x) 
dgdx = x*(6 .  *c3  + 1 2 .  *x) 
xl=x-g/dgdx 

if  ((abs(g).gt.0.0000001).or.(abs(xl-x).gt.0.0000001)) 
+  goto  10 

C-  Now  apply  Newton's  method  to  find  the  roots 
f=rcfunc(xl,c0,cl,c3) 
if  (f.gt.0.)  then 
rc=-998 
else 

do  1  j  =  l,2 
if  (j.eq.l)  then 
mult=-l 


20 


1 


else 

mult=l 
end  if 

x2=xl*l.+0.1*mult 


x=x2 

f0=rcfimc(x,c0,cl  ,c3) 
fl=cl+x*x*(3.*c3+4.*x) 
dxl=-fO/fl 
x2=dxl+x 

if  (abs(l.-(x2/x)).gt.0.0001)  goto  20 
if  (mult. It. 0)  rc=x2 
if  (mult.gt.0)  re=x2 
continue 
end  if 
end  if 
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return 

end 


Real*8  Function  rcfunc(x,c0,cl,c3) 

implicit  real*8  (a-h,o-z) 

z=c0+x*(cl+x*x*(c3+x)) 

rcfunc=z 

return 

end 


c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

subroutine  regrid(ia,aa,gridsize,stepsize) 
implicit  real*8  (a-h,o-z) 
dimension  ia(3,24),aa(3,24) 
do  1901  iii=  1,3 
do  1901  iiii=l,24 

aa(iii,iiii)=dfloat(ia(iii,iiii))*gridsize 
1901  continue 

stepsize = gridsize*dsqrt(3 .  dO) 
return 
end 
c 
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Section  D5.  Fortran  source  code  for  Monte  Carlo  simulation  which  utilizes  radially  symmetry  of  problem  to 
speed  computation. 

C-  This  program  simulates  bubble  growth  beginning  with  a  bubble  of 
C-  critical  radius  as  defined  by  Tikuisis.  The  bubble  is  represented  as 
C-  a  region  of  high  solubility  surrounded  by  a  "depleted  region"  from 
C-  which  the  inert  gas  is  extracted.  These  two  regions  are  placed  in 
C-  a  third  region  which  represents  the  surrounding  tissue, 
c  this  program  also  uses  1  dimension,  the  radial  one,  to  keep 
c  track  of  the  particle,  should  speed  things  up  considerably. 

implicit  real*8  (a-h,o-z) 
common/com2/pi,  i  1 
common/comroot/probroot,firstime 
common/distarr/distance( 1 0000) ,  sd(500) 
common  /com5/twopi,onethird,vcoeff 
real*8  N2,N2dep,k,kH,molezone(500) 

Dimension  xx(3),sumcumm(0:500) 
dimension  endzone(500),fraczone(500),conczone(500) 
dimension  cummsum(0:500),rzmgb(500),rzrnge(500) 
dimension  rzb(0 : 500) ,  rz2b(0 : 500) ,  rz3b(0 : 500) ,  vzb(500) 
dimension  rze(0 : 500) ,  rz2e(0 : 500) ,  rz3e(0 : 500) , vze(500) 

C-  Define  constants 
zed=0.d0 

Pi = ■ 4 .  dO*datan(  1 .  dO) 
onethird  =  1 .  dO/3 .  dO 
vcoeff = 4 .  dO*pi*onethird 
twopi=2.d0*pi 

c  change  12/16/92  to  include  oxygen  and  co2  as  additional  gas 
c  components  of  the  bubble,  using  values  from  van  liew,  pg.  336 
c 

c  ph2o  =  6.2  kPa  =  46.52  mmHg 

c  po2  =5.3  kPa  =  39.76  mmHg 

c  pco2  =5.9  kPa  =  44.26  mmHg 

c  pco21ung  =  5.3  kPa  =  39.76  mmHg 
c  pvap=  46.52d0  +  39.76d0  +  44.26d0 

c  read  in  from  input  file  below  1/18/93 
c  ph2o=46.52d0 

c  po2=39.76d0 

c  pco2=44.26d0 

c  pco21ung=39.76d0 

c  pvap=ph2o+po2+pco2 

C-  vapor  pressure  of  water  in  bubble  @  37  C  in  mmHg 
C1125 =0.055346 
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C11=0. 055140 

C-  pure  solvent  cone  of  water  @  37  C  in  moles/cc 
k=62358.0 

C-  Boltzman's  gas  constant  in  mmHG-cc/mole-K 
C-  based  on  PV=kNT  and  STP  conditions 
gasR= 62358.0 

C-  Universal  gas  constant  in  mmHg-cc/mole-K 
C-  based  on  published  value  of  0.08205  L-atm/mole-K 
T=37. 0+273.1 
C-  temperature  in  Kelvin 

C-  Input  initial  conditions  and  control  variables. 

C-  ngrow= number  of  bubble  growth  cycles, 

C-  ntrip= number  of  independent  trips  thru  tissue,  nstep= number  of  steps  per 
C-  trip,  Surten= surface  tension  of  bubble  (dynes/cm),  PCbub = partition 
C-  coefficient  of  bubble  (1/ostwald  coefficient),  D= diffusion  coefficient 
C-  (cm  sq/sec),  Voldep0= initial  volume  of  the  depleted  region  (cu  cc), 

C-  Voltiss  =  volume  of  cubical  tissue  region 

C-  P01= initial  tissue  pressure  (mmHg),  Pstep= decompression  step  (mmHg) 
Rewind(l) 

Read(l,900)ngrow,ntrip,nstep,surten,Pcbub,D,gridsize,Voldep0, 

+  Voltiss, P01,Pstep,bubfact 

900  Format(I12/I12/I12/gl5.7/gl5.7/gl5.7/gl5.7/gl5.7/ 

+  g15.7/gl5.7/gl5.7/gl5.7) 

read(l ,  1900)iprofile,nmicro,washoutt,pctn2 

1900  format(il2/il2/gl5.7/gl5.7) 

read(l ,  1 90  l)s  witch,  volref,ffacsd,ph2o,po2,pco2,pco21ung 

1901  format(gl5.7) 
close(l) 

pvap =ph2o +po2 +pco2 

C-  Input  the  array  of  diffusion  distances  as  a  function  of  probability 
C-  given  in  terms  of  the  standard  deviation 

Open(unit=  1  ,file  =  'distprob.dat'  ,type  =  'old') 

Read( 1 , 899)(distance(i) ,  i  =  1 , 10000) 

899  Format(fl7.2) 

Close(l) 

vol = voltiss + voldepO 

C-  adjust  units  to  mmHG-cm  knowing  760  mmHg=  1.013X10*6  dyn/cnU2 
Surten = Surten/1 332 . 8 

C-  Calculate  Henry's  coefficient  from  PCbub 

KH = 760 .  /((((gasR*T)/((  1 .  /c  1 1)  *(  1 .  /PCbub)  *760 . )) + 1 )  *  *- 1 . ) 

C-  Calculate  stepsize 
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stepsize = gridsize  *dsqrt(3 .  dO) 

C-  Calculate  bubble  growth  cycle  time  and  tissue  washout  time 
Dtime = stepsize  *steps  ize/(6.  *D) 
eyelet = nstep  *Dtime 
do  333  kk=l, nstep 
xtime=kk*dtime 
sd(kk) = dsqrt(2  .dO*d*xtime) 

333  continue 

poofprob = 1  .dO-dexp(-cyclet/washoutt) 
if(washoutt .  gt .  86400 .  d0)poofprob = 0 .  dO 
C-  Calculate  the  maximum  distance  a  particle  can  travel  by  diffusion 
onesd = dsqrt(2  .dO*d*cyclet) 
diffdist = 5 .  d0*onesd 
c 

c  figure  out  the  width  of  the  diffusion  shells,  switch  =  0  is  for 
c  a  dimensionless  run,  with  reference  volume  volref.  switch  <  >  0  is 
c  an  absolute  run,  with  parameters  unsealed. 

c 

if  (switch .  eq .  0 .  d0)then 
width = onesd*fracsd*(vol/volref)*  *onethird 
if(width .  gt .  onesd)  width = one  sd 
else 

width = onesd*ff  acsd 
end  if 

if(iprofile .  eq .  0)then 
C-  Write  out  the  starting  conditions 

Open(unit=2,file  =  'for002.dat' ,type=  'new') 

Write(2,901)  Ngrow,Ntrip,Nstep,Surteh*1332.8,PCbub,D, 

+  gridsize ,  stepsize ,  P01 ,  Pstep ,  VoldepO ,  voltiss 

Close(2) 

901  Format(3I9/5G15.4/4G15.4) 

end  if 

C-  Begin  a  bubble  growth  cycle 
time=0. 

P1=P01 

pln2 = (p01-ph2o-pco21ung)*pctn2 

pfn2 = (p01-pstep-ph2o-pco21ung)*pctn2 

Voldep= VoldepO 

il  =  182361  +2*secnds(0.0) 

lstrt=l 

if (iprofile .  ne .  0)then 

c  for  a  nonuniform  initial  concentration  profile,  put  code  in  here 

open(unit=7,file=  'bubmod.dmp',type  =  'old',form=  'unformatted') 
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rewind(7) 

read(7)lstrt,il, nzonebeg, time, volbub,calcpc,rbub,rdep,rbubO 
read(7)re ,  rc ,  currmole 

read(7)(ffaczone(i),molezone(i),conczone(i),i  =  1  ,nzonebeg) 

read(7)(rzb(i) ,  rz2b(i) ,  rz3b(i) ,  vzb(i) ,  rzmgb(i) ,  i  =  1 ,  nzonebeg) 

read(7)(rze(i),rz2e(i),rz3e(i),vze(i),rzmge(i),i= 1  ,nzoneend) 

close(7) 

lstrt=lstrt+ 1 

cs=cll*pln2/kh 

N2=Cs*Vol 

N2dep = Cs*Voldep 

orgmoles=N2 

If  (VoldepO.gt.vol)  orgmoles=n2dep 

Pl=Pl-Pstep 

csequil = c  1 1  *pfn2/kh 

rzb(0)=0.d0 

rze(0)=0.d0 

bubmoles =molezone(l) 

open(unit=3,file=  'for003.dat',  type= 'new') 

write(3 , 3000)0 ,  time ,  rbub ,  volbub ,  calcpc , 

+  (fraczone(i),rzb(i),conczone(i),i=  1 , nzonebeg) 

close(3) 

depmoles = currmole-bubmoles 
bubff  ac = bubmoles/currmole 
depff ac = depmoles/currmole 
end  if 

Do  10  L=lstrt,ngrow 

C-  Calculate  number  of  moles  of  inert  gas  (N2)  at  start  of  the  cycle 
C-  in  the  entire  tissue  and  in  the  depleted  region  alone  (N2dep). 

If((L.eq. l).and.(iprofile.eq.O))  then 
5432  continue 

cs=cll*pln2/kh 
N2=Cs*Vol 
N2dep = Cs*Voldep 
orgmoles=N2 

If  (VoldepO.gt.vol)  orgmoles=n2dep 
currmole = orgmoles 
Pl=Pl-Pstep 
csequil = c  1 l*pfn2/kh 
C-  Calculate  bubble  dimensions 
C-  First,  calculate  the  critical  and  equilibrium  radii 
aO = 2 .  *surten*Voldep/Pl 
Cs=Cll*Pl/KH 
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al  =  Voldep*(l  .-Pvap/Pl)-N2dep/Cs 

a3 = 8 .  *pi*surten/(3 .  *Cs  *k*T) 

a4 = 4 ,  *pi  *P1  *(  1 .  -Pvap/Pl)/ (3 .  *Cs  *k*T) 

c0=a0/a4 

cl=al/a4 

c3=a3/a4 

Call  roots(cO,cl,c3,rc,re) 
if  (rc.lt.O.)  goto  999 
C-  Assign  radius  to  test  bubble 

rbub = re*bubfact + rc  *(  1 .  dO-bubfact) 
rbub =dmaxl(rbub,(l  .5d0*stepsize)) 
rbubO=rbub 

VolBub =4.  *pi*rbub*rbub*rbub/3 . 
rdep = (3 .  *(Voldep + Volbub)/(4*pi))**(l  ./3 .) 
c  check  here  to  make  sure  the  bubble  fits  inside  the  cube,  if  not,  spread 
c  the  excess  volume  onto  the  depleted  region,  and  make  it  larger 
rzb(0)=0.d0 
rze(0)=0.d0 

call  fastgenz(rbub, rdep, width, nzonebeg,rzb,rz2b,rz3b,vzb,rzrngb) 
Bubmoles = Volbub*(Pl-Pvap + 2 .  *surten/rbub)/(gasR*T) 

Depmoles = N2Dep-Bubmoles 
Tismoles = N2-N2dep 
Bubffac =Bubmoles/orgmoles 
Depfrac = Depmoles/orgmoles 
Tisffac = Tismoles/orgmoles 
if(iprofile.eq.O)then 
ff  aczone(  1 ) = bubfrac 
molezone(  1 ) = bubmoles 
conczone(l) =bubmoles/volbub 
do  1100  izone=2,nzonebeg 

fraczone(izone) = depfr  ac*vzb(izone)/voldep 
molezone(izone) = fraczone(izone)*orgmoles 
conczone(izone) = depmoles/voldep 
1 100  continue 
end  if 
End  if 

call  fastgenz(rbub, rdep, width, nzoneend,rze,rz2e,rz3e,vze,rzmge) 

C-  Begin  the  random  walk  thru  the  module 
C-  Calculate  the  probability  of  staying  in  the  bubble 
rbpw = rbub + width 
stepsiz2 = stepsize*stepsize 
rbub2 = rbub  *rbub 
rbub3  =rbub2*rbub 
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rdep2=rdep*rdep 

pf  =  stepsize  *(  12  .d0*rbub2-stepsiz2)/(  1 6 .  d0*(rbub3)  *pcbub) 
ddprbub = diffdist + rbub 
firstime=zed 
EndInbub=0. 

EndIndep=0. 

EndIntis=0. 

do  1005  izone=l,nzoneend 
endzone(izone) = O.dO 
1005  continue 

sumcumm(O) = zed 
cummsum(O) = 0 .  dO 
do  1004  izone  =  l,nzonebeg 

cummsum(izone) = cummsum(izone-l) + fraczone(izone) 

1004  continue 

if(cummsum(nzonebeg).gt.zed)then 
if(voltis .  le .  zed)then 
do  1904  izone=0,nzonebeg 

cummsum(izone) = cummsum(izone)/cummsum(nzonebeg) 
sumcumm(nzonebeg-izone) = 1  .dO-cummsum(izone) 

1904  continue 
end  if 
end  if 

C-  Calculate  the  number  of  particles  to  be  place 
Do  20  n=l,Ntrip 

C-  Place  a  particle  randomly  in  the  module  with  probability  weighted 
C-  by  mole  fraction. 

prob=ran(il) 

c 

c  new  code  to  place  particles  in  the  various  zones 
do  1001  izone  =  l,nzonebeg 
if(prob .  It.  sumcumm(izone))then 
call  place(r,rz3b(nzonebeg-izone),rzmgb(nzonebeg-izone  + 1)) 
r2=r*r 
rp=r 
rp2=r2 
goto  1002 
end  if 

1001  continue 

1002  continue 

C-  If  the  particle  is  placed  so  far  from  the  bubble  that  it  cannot 
C-  reach  it  in  nstep  steps  and  the  module  only  consists  of  the 
C-  bubble  and  a  depleted  region,  the  particle  will  stay  in  the  depleted 
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C-  region. 

If  (r.gt.ddprbub)  then 
Istep=nstep 
call  getsize(a,a2,istep) 
call  faststep(a,a2,r,r2,rp,rp2,cosphi) 
goto  31 

Else  If  (r.gt.rbpw)  then 

C-  The  particle  can  reach  the  bubble  by  diffusion 
Istep = int((r-rbub)/stepsize) 
call  getsize(a,a2,istep) 
call  faststep(a,a2,r,r2,rp,rp2,cosphi) 

Else 

C-  The  particle  is  in  the  bubble 
Istep =0 
End  if 

C-  Carry  out  random  walk  for  NSTEP  steps  of  "stepsize" 

a=stepsize 

a2=stepsiz2 

Do  30  JJ=Istep+l,nstep 
lnbub=0 

C-  Check  to  see  if  particle  is  in  the  bubble  by  comparing  the 
C-  particle  distance  from  origin  with  bubbble  radius 

if(rp.lt.rbub)inbub = 1 

C-  If  the  particle  is  in  the  bubble,  leave  it 
C-  based  on  the  probability  factor. 

If  (Inbub.eq.l)  then 
if  (ran(il).lt.pf)then 
call  fastout(r,a,rbub) 
rp=r 

rp2=rp*rp 
end  if 
else 

C-  Take  a  step 
r=rp 
r2=rp2 

call  faststep(a,a2,r,r2,rp,rp2,cosphi) 
if(rp.gt.rdep)then 
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call  fastrefl(a ,  r ,  r2  ,rp ,  rp2 ,  rdep  ,rdep2 ,  cosphi) 
r=rp 
r2=rp2 
end  if 
End  if 

30  Continue 

31  Continue 
if(rp.gt.rdep)then 

call  fastrefl(a,r,r2,rp,rp2, rdep, rdep2, cosphi) 

r=rp 

r2=rp2 

if(r.gt.rdep)then 

c  if  the  particle  is  still  outside  the  depletion  region  after  reflection, 
c  assume  that  further  reflections  will  leave  it  in  the  outermost  zone, 
call  place(r ,  rz3b(nzonebeg- 1 )  ,rzmgb(nzonebeg)) 

r2__r*r 

end  if 
else 
r=rp 
r2=rp2 
end  if 

C-  Keep  track  of  the  region  in  which  the  particle  ends 
do  1010  izone=l,nzoneend 
if(r .  It.  rze(izone))then 

endzone(izone) =endzone(izone)  + 1  .d0 
goto  1011 
end  if 

1010  continue 

101 1  continue 

If  (r.le.rbub)  then 
Endlnbub = Endlnbub  + 1 . 
else  If  O',  le.  rdep)  then 
Endlndep = Endlndep + 1 . 
else 

Endlntis = Endlntis  + 1 . 
end  if 

20  Continue 

if((voltiss .  le .  zed) .  and .  (endintis .  gt .  zed))then 
open(unit=2,file  =  'for002.dat'  ,type  =  'old'  .access  =  'append') 
write(2,8888)endintis  ,1 

8888  format('  accounting  error... ',fl0.2,' particles  in  tissue'/ 

+  '  on  the  ',i8,'th  cycle') 


117 


stop 
end  if 

C-  Calculate  the  gas  content  in  each  region 
If  (L.eq.l)  then 

open(unit=2,file  =  'for002.dat' , type  =  'old' .access  =  'append') 
write(2,904)time,rbub,volbub,voldep, 

+  (vol-voldep) 

write(2,904)cyclet,poofprob,washoutt 
Write  (2,905)  rc,re,rbubO 
close(2) 
end  if 

dfnt = dfloat(ntrip) 
dfntO=dfht 
actual =endindep 
if(poofprob.gt.0.d0)then 
const = dfht*csequil/currmole 
expected = const*  voldep 

actual = endindep-(endindep-expected)  *poofprob 
dfnt = actual + endinbub 
1025  continue 
1015  continue 

if(voltiss .  gt .  zed)dfnt = dfnt + endintis 
currmole = currmole*dfht/dfhtO 
end  if 

do  1020  izone=l,nzoneend 

fraczone(izone) = endzone(izone)/dfht*actual/endindep 
if(izone.eq.  l)ffaczone(izone) =endzone(izone)/dfht 
molezone(izone) = fraczone(izone)*currmole 
conczone(izone) = molezone(izone)/vze(izone) 

1020  continue 

calcpc = fraczone(l)*vze(2)/(ffaczone(2)*vze(l)) 

Bubfrac = EndlnBub/dfnt 
depfrac = actual/dfnt 
Tisffac = Endintis/dfnt 
bubmoles = currmole  *bubffac 
depmoles = currmole  *depfrac 
tismoles = currmole  *tisffac 

C-  If  Bubfrac  is  very  small  it  will  soon  vanish  so  trap  it 
If  (bubfrac. lt.0. 0000000001)  goto  999 

C-  Calculate  the  new  volume  of  each  region 

C-  the  volume  of  the  depleted  region  plus  the  volume  of 
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C-  the  supersaturated  tissue  region  is  a  constant,  so  bubble 
C-  growth  only  adds  to  max  and  min  dimensions,  but  not  volumes 
A3 = (4.  *pi*(Pl-Pvap))/3 . 

A2 = (8 .  *pi*surten)/3 . 
aO = -molezone(l)*gasR*T 
CO = AO/  A3 
C2=A2/A3 

Call  Roots2(c0,c2,rbub) 
if(rbub.lt.stepsize)goto  999 
VolBub = (4*pi*(rbub)**3  .)/3 . 
tis  VNew = T  ismoles/(C  11*P01/KH) 

Voldep = Vol-tisvnew 

If  (VoldepO.ge.vol)  Voldep =voldepO-tisvnew 
rdep = ((voldep + volbub)/vcoeff)**one  third 
call  fastexpnd(rbub,rdep,nzoneend,rzb,rz2b,rz3b,rzrnge,rzrngb) 
nzonebeg = nzoneend 

C-  Output  information  for  1  bubble  growth  cycle 
time = time + eyelet 

open(unit=2,file  =  'for002.dat'  ,type  =  'old' , access  =  'append') 

Write(2,904)  time,rbub,calcpc,currmole 

close(2) 

904  Format(5G15.7) 

3000  formate  cycle  #  ',i5,4gl2.4/50(3gl5.6/» 

open(unit=7,file=  'bubmod.dmp'  ,type  =  'old'  ,form=  'unformatted') 
rewind(7) 

write(7)L,  il ,  nzonebeg,  time ,  volbub  ,calcpc  ,rbub ,  rdep  ,rbub0 
write(7)re,rc,currmole 

write(7)(ffaczone(i),molezone(i),conczone(i),i= 1  .nzonebeg) 
write(7)(rzb(i),rz2b(i),rz3b(i),vzb(i),rzrngb(i),i  =  1  .nzonebeg) 
write(7)(rze(i),rz2e(i),rz3e(i),vze(i),rzmge(i),i = 1  .nzoneend) 
close(7) 

10  Continue 
999  Continue 

905  Format(3G15.4) 


Stop 

End 


c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

subroutine  fastout(r,a,rb) 
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implicit  real*8  (a-h,o-z) 
common  /com2/pi,il 
common/comroot/probroot,firstime 
probroot=ran(il) 
s=rtsafe(0.d0,l.d0,l.d-10,rb,a)*a 
r=s+rb 
return 
end 
c 
c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  subroutines  to  generate  c.d.f.  for  particles  stepping  out  of  a  bubble, 
c  some  parameters  are  supplied  to  the  subroutine  via  the  common  block 
c  comroot. 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

FUNCTION  RTSAFE(Xl,X2,XACC,b,a) 
implicit  real*8  (a-h,o-z) 

PARAMETER  (MAXIT=100) 

CALL  FUNCD(X1  ,FL,DF,b,a) 

CALL  FUNCD(X2,FH,DF,b,a) 

IF(FL*FH.GE.O.)  PAUSE  'root  must  be  bracketed' 

IF(FL.LT.O.)THEN 

XL=X1 

XH=X2 

ELSE 

XH=X1 

XL=X2 

SWAP=FL 

FL=FH 

FH=SWAP 

ENDIF 

RTSAFE=  .5*(X1  +X2) 

DXOLD=ABS(X2-Xl) 

DX=DXOLD 

CALL  FUNCD(RTSAFE,F,DF,b,a) 

DO  11  J  =  1  ,MAXIT 

IF(((RTSAFE-XH)*DF-F)*((RTSAFE-XL)*DF-F).GE.O. 

*  .OR.  ABS(2.  *F).GT.  ABS(DXOLD*DF) )  THEN 
DXOLD=DX 
DX=0.5*(XH-XL) 

RTS  AFE = XL + DX 
IF(XL.  EQ  .RTS  AFE)RETURN 
ELSE 
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DXOLD=DX 
DX=F/DF 
TEMP = RTS  AFE 
RTS  AFE = RTS  AFE-DX 
IF(TEMP.  EQ .  RTS  AFE)RETURN 
ENDIF 

IF(ABS(DX).LT.XACC)  RETURN 
CALL  FUNCD(RTSAFE,F,DF,b,a) 

IF(F.LT.O.)  THEN 
XL = RTS  AFE 
FL=F 
ELSE 

XH=RTSAFE 

FH=F 

ENDIF 

11  CONTINUE 

PAUSE  'RTS AFE  exceeding  maximum  iterations' 

RETURN 

END 

c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

subroutine  funcd(x,f,d,b,a) 
implicit  real*8  (a-h,o-z) 
common/comroot/probroot,firstime 
if(firstime.eq.0.dO)then 
roa=b/a 
roa2=roa*roa 
f  1  =  (1 2 .  d0*roa-24 .  d0*roa2) 
f2 = (6.d0-24.d0*roa + 12.d0*roa2) 
f3=(12.dO*roa-8.dO) 
f4=3.d0 
dO=fl 
dl=2.d0*f2 
d2=3.d0*f3 
d3=4.d0*f4 
firstime=l.dO 
end  if 

f0= (12.d0*roa2-l  ,dO)*probroot 

f=((((f4*x+f3)*x+f2)*x+fl)*x+fD) 

d=(((d3*x+d2)*x+dl)*x+d0) 

return 

end 
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cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

subroutine  fastexpnd(rb,rd,nz,rz,rz2,rz3 ,rzre,rzrb) 
implicit  real*8  (a-h,o-z) 

dimension  rz(0 : 500) ,  rz2(0 : 500) ,  rz3 (0 : 500) ,  rzre(500) ,  rzrb(500) 
common/com5/twopi,onethird,coeff 
data  nzmax,eps/500, 1  .d-5/ 
rz(0)=0.d0 
rz2(0)=0.d0 
rz3(0)=0.d0 
rz(l)=rb 
rz2(l)=rb*rb 
rz3(l)=rb*rz2(l) 
rzrb(l)=rzre(l) 
do  10  i=2,nz 
rz3(i) =rz3(i-l) +rzre(i) 
rz(i) =rz3(i)**onethird 
rz2(i)=rz(i)*rz(i) 
rzrb(i)=rzre(i) 

10  continue 

do  50  i=nz+l,nzmax 
rz(i)=0.d0 
rz2(i)=0.d0 
rz3(i)=0.d0 
rzrb(i)=rzre(i) 

50  continue 

if(dabs(rz(nz)-rd) .  gt .  eps)then 
write(3 , 100)nz,rz(nz),rd 

100  format(’  expansion  error:  rz(',i3,')  =  ',gl5.6, 

+  '  rdep  =  \gl5.6) 

stop 
end  if 
return 
end 
c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 


c  variables: 

a 

=  stepsize 

c 

rl 

=  initial  radial  position 

c 

rl2 

=  rl*rl 

c 

r2 

=  radial  distance  >  rd  (do  not  use) 

c 

r22 

=  r2*r2  (do  not  use) 

c 

r3 

=  final  radial  position 

c 

r32 

=  r3*r3 
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c 

rd 

=  depletion  radius 

c 

rd2 

=  rd*rd 

c 

cp 

=  cos(phi),  from  faststep 

c 

St 

=  sin(theta) 

c 

st2 

—  st*st 

c 

ct 

=  cos(theta) 

c 

c  NOTE:  even  after  reflection  from  the  outer  surface,  the  particle 
c  may  still  be  outside  the  depletion  region,  have  to  check  for 

c  this  in  the  main  program, 

c 

subroutine  fastrefl(a,rl ,rl2,r3,r32,rd,rd2,cp) 

implicit  real*8  (a-h,o-z) 

common/com2/pi,  i  1 

common/com5/twopi,onethird,coeff 

cp2=cp*cp 

sp2=l.d0-cp2 

if(sp2 .  It .  0 .  d0)sp2 = 0 .  dO 

dl  =rl*cp+dsqrt(rd2-rl2*sp2) 

St2=rl2*sp2/rd2 

if(st2.ge.l.d0)then 

ct=0.d0 

else 

ct=dsqrt(l.d0-st2) 
end  if 
d2=a-dl 
d22=d2*d2 

r32 = rd2 + d22-2  ,d0*rd*d2*ct 
r3  =dsqrt(r32) 
return 
end 
c 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCGCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

subroutine  faststep(a ,  a2 ,  r ,  r2 ,  rp ,  rp2 ,  cosphi) 
implicit  real*8  (a-h,o-z) 
common/com2/pi,il 
common/com5/twopi,onethird,coeff 
cosphi = 1  .d0-2.d0*ran(il) 
rp2 = a2 + r2-2  .dO*a*r*cosphi 
rp=dsqrt(rp2) 
return 
end 
c 
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cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


subroutine  place(r,rz3,rzmg) 
implicit  real  *8  (a-h,o-z) 
common/com2/pi,il 
common/com5/twopi,onethird,coeff 
r = (rz3  +  rzmg  *ran(i  1 ))  *  *onethird 
return 
end 
c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

subroutine  fastgenz(rb,rd,w,nz,rz,rz2,rz3,vz,rzmg) 
implicit  real*8  (a-h,o-z) 

dimension  rz(0:500),rz2(0:500),rz3(0:500),vz(500),rzrng(500) 

common/com5/twopi,onethird,coeff 

data  nzmax/500/ 

nz=  1 

rz(0)=0.d0 
rz(l)=rb 
do  10  i=2,nzmax 
nz=nz+l 
rz(i)=rz(i-l)+w 
if(rz(i).ge.rd)then 
rz(i)=rd 
goto  20 
end  if 
10  continue 

20  continue 

rz2(0)=0.d0 
rz3(0)=0.d0 
do  40  i=l,nz 
rz2(i)=rz(i)*rz(i) 
rz3(i)=rz2(i)*rz(i) 
rzmg(i) = rz3(i)-rz3(i-l) 
vz(i) = coeff'rzmg(i) 

40  continue 

do  50  i=nz+l,nzmax 
rz2(i)=0.d0 
rz3(i)=0.d0 
rzmg(i)=0.d0 
vz(i)=0.d0 
50  continue 

return 
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end 


c 

Subroutine  Roots(c0,cl,c3,rc,re) 
implicit  real*8  (a-h,o-z) 

C-  This  subroutine  calculates  the  roots  of  the  Tikuisis  polynomial. 
C-  First,  find  the  minimum  because  1  root  is  below  it  and  the  other 
C-  root  is  above  it. 

100  formate  #  ’45) 

If  (cl.ge.O.)  then 
rc=-999. 

else  ' 

xl=-c0/cl 
10  x=xl 

g = c  1  +  x*x*(3 .  *c3  +  4 .  *x) 
dgdx = x*(6 .  *c3  + 1 2 .  *x) 
xl=x-g/dgdx 

if  ((abs(g).gt.0. 0000001). or. (abs(xl-x).gt. 0.0000001)) 

+  goto  10 

C-  Now  apply  Newton's  method  to  find  the  roots 
f =rcfixnc(xl  ,c0,cl  ,c3) 
if  (f.gt.0.)  then 
rc=-998 
else 

do  1  j  =  l,2 
if  (j.eq.l)  then 
mult=-l 
else 

mult=l 
end  if 

x2=xl*l.+0.1*mult 
20  x=x2 

fO =rcfunc(x,c0,cl  ,c3) 
f  1  =  c  1  +  x*x*(3 .  *c3  +4 .  *x) 
dxl=-f0/fl 
x2=dxl+x 

if  (abs(l.-(x2/x)).gt.0.0001)  goto  20 
if  (mult.  It.  0)  rc=x2 
if  (mult.gt.0)  re=x2 
1  continue 

end  if 
end  if 
return 
end 
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Real*8  Function  rcfunc(x,c0,cl,c3) 

implicit  real*8  (a-h,o-z) 

z = cO + x*(c  1  +  x*x*(c3  4-  x)) 

rcfiinc=z 

return 

end 


£********************************  ***  ************************************ 
Subroutine  Roots2(c0,c2,rbub) 
implicit  real*8  (a-h,o-z) 

C-  This  subroutine  calculates  the  roots  of  the  polynomial  describing 

C-  bubble  radius. 

C-  Apply  Newton's  method  to  find  the  roots 
x2 = rbub 

20  x=x2 

fD =rbubfunc(x,c0,c2) 
fl=x*(2.*c2+3.*x) 
dxl=-fO/fl 
x2=dxl+x 

if  (abs(l.-(x2/x)).gt.0.0001)  goto  20 
rbub=x2 
return 
end 


Real*8  Function  rbubfunc(x,c0,c2) 
implicit  real*8  (a-h,o-z) 
z=c0+x*x*(c2+x) 
rbubfimc=z 
return 
end 
c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

c  subroutine  to  take  find  the  size  of  a  single  large  diffusion 
c  step  for  use  in  subroutine  faststep 

c 

c  subroutine  getsize(a,a2,diffcoef,t) 

subroutine  getsize(a,a2,isd) 
implicit  real*8  (a-h,o-z) 
common  /com2/pi,il 
common/distarr/distance( 1 0000) ,  sd(500) 
c 

c  the  distance  stepped  as  a  function  of  the  probability,  given  in  terms 
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c  of  the  standard  deviation,  is  located  infile  distprob.dat 
c 

index = aint(  1 .  dO + 1 0000 .  dO*ran(i  1 )) 
a = sd(isd)  *distance(index) 
a2=a*a 
return 
end 
c 
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